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ABSTRACT 

Geometric  Algebra  is  the  most  appropriate  unifying  mathematical  language 
to  describe  diverse  problems  in  mathematics,  physics,  engineering  and  com¬ 
puter  science.  In  combination  with  Projective  Geometry  it  provides  an  effi¬ 
cient  framework  for  computer  vision  and  robotics,  where  image  processing  and 
recognition  play  the  central  role.  This  document  addresses  a  gentle  introduc¬ 
tion  to  Geometric  Algebra  followed  by  its  implementation  in  MATLAB.  The 
developed  fully  vectorized  code  is  thoroughly  tested.  Several  applications  are 
presented  in  the  form  of  unit  tests,  amongst  which  are  some  basic  algorithms 
for  the  reconstruction  of  a  three-dimensional  structure  from  two-dimensional 
images. 
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1  Introduction 

It  has  been  recognized  that  Geometric  Algebra  is  the  most  appropriate  unifying  mathe¬ 
matical  language  to  describe  diverse  problems  in  mathematics,  physics,  engineering  and 
computer  science  [Hestenes  &  Sobczyk  1987,  Corrochano  <fc  Sobczyk  2001,  Perwass  2009]. 
In  particular,  in  combination  with  Projective  Geometry  it  provides  an  efficient  framework 
for  computer  vision  and  robotics,  where  image  processing  and  recognition  play  the  central 
role. 

There  are  several  algebraic  structures  suitable  for  the  formalization  of  geometric  oper¬ 
ations,  such  as  Gibbs  vector  algebra,  Grassmann  algebra,  Grassmann-Cayley  algebra  and 
Clifford  algebra.  The  classical  vector  algebra  operates  on  vectors  but  has  limitations  when 
dealing  with  more  complex  geometric  constructs.  Grassmann  algebra  is  equipped  with 
an  exterior  product,  also  called  a  progressive  outer  product,  which  is  useful  for  building 
higher  dimensional  objects  from  lower  dimensional  objects,  for  example,  when  creating 
a  line  joining  two  points  or  a  plane  passing  through  three  points  in  a  projective  space. 
Grassmann-Cayley  algebra  is  additionally  equipped  with  a  regressive  outer  product  which 
formalizes  the  operation  of  intersection  of  linear  objects,  such  as  the  creation  of  a  line  of 
intersection  between  two  planes.  Finally,  Clifford  algebra  additionally  contains  an  inner 
product  which  is  used  to  describe  parallelism  and  orthogonality. 

In  this  document  we  confine  ourselves  to  Clifford  algebra  generated  by  a  vector  space 
equipped  with  a  positive  definite  scalar  product.  Note  that  the  Grassmann  and  Grassmann- 
Cayley  algebras  can  be  factored  out  from  Clifford  algebra,  and  hence  they  are  part  of  it. 
In  our  context  these  algebraic  structures  are  identified  with  Geometric  Algebra. 

We  provide  a  gentle  introduction  to  Geometric  Algebra  and  Projective  Geometry  to 
fix  definitions  and  notations,  and  recall  basic  relations  between  algebraic  and  geometric 
operations.  The  developed  MATLAB  code  for  a  range  of  Clifford  algebras,  based  on  the 
algebra  multiplication  tables  given  in  Appendix  A,  is  presented  in  Appendix  B.  The  fully 
vectorized  code  is  thoroughly  tested  against  the  basic  properties  of  the  Clifford  algebra 
operations  and  derived  identities  by  a  suite  of  unit  tests  given  in  Appendix  C. 

When  combined  with  projective  spaces,  the  real  power  of  Clifford  algebra  emerges.  In 
particular,  the  effectiveness  of  Geometric  Algebra  manifests  itself  by  an  elegant  and  concise 
form  of  equations  for  intricate  geometric  constructs.  We  present  another  suite  of  unit  tests 
in  Appendix  D  which  demonstrates  this  power.  We  apply  this  technique  to  the  derivation 
of  basic  algorithms  for  the  projection  of  a  three-dimensional  object  to  two-dimensional 
images,  followed  by  its  reconstruction  from  a  sequence  of  images. 


2  Geometric  algebra 

Consider  an  n-dimensional  vector  space  Mn  over  the  field  of  real  numbers  M,  equipped 
with  a  positive  definite  scalar  product  denoted  by  the  centred  dot 

Mn  x  Mn  3  (a,b)  a  •  b  €M..  (1) 

Choose  an  orthonormal  basis  {ei,...,en}  characterized  by  the  identities  e,  •  ej  =  8ij 
where  5ij  is  the  Kronecker  symbol.  The  Clifford  algebra  C  (Mn)  is  defined  by  the  following 
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multiplication  rules  on  the  basis  elements 

ci  Cj  A  Cj  e*  —  26ij  (2) 

in  combination  with  the  additional  axioms  of  associativity  and  distributivity  [Lang  2002, 
Gar  ling  2011,  Hestenes  &  Sobczyk  1987]. 

We  denote  the  multiplication  operation  in  the  Clifford  algebra  by  juxtaposition.  This 
product  is  called  the  geometric  or  Clifford  product  to  distinguish  it  from  other  products 
to  be  introduced.  The  multiplication  rules  on  the  basis  vectors  can  be  written  as 

Cj  Cj  —  •  Cj  +  Cj  A  c,j  ( .3  ) 

where  the  inner  (scalar)  and  outer  (wedge)  products  are  expressed  for  any  vectors  a,  b  in 
terms  of  the  geometric  product 

a  ■  b  =  -  (a  b  +  b  a)  ,  a  A  b  =  -  (a  b  —  b  a)  .  (4) 

For  the  basis  orthonormal  vectors  we  have  ef  =  1  and  Cj  Cj  =  Cj  A  Cj  =  —  Cj  A  e*  =  —  ej  Cj 
provided  that  i  j -  j.  The  elements  of  C  (Rn)  are  called  multivectors,  hypercomplex  numbers 
or  Clifford  numbers. 

The  basis  multivectors  of  the  Clifford  algebra  are  constructed  by  the  collection  {e/} 
where  the  compound  index  I  is  a  subset  of  the  set  {1,  2, . . . ,  n}.  For  convenience  we  adopt 
the  notation  e@  =  1  (when  I  =  0)  and  set 

6/  22  ...ik  ~~  ©l  &i2  •  *  •  Cj^  for  I  ^2?  ■  •*?  (b) 

where  we  will  always  assume  that  1  <  i\  <  12  <  •  •  •  <  ijt  <  n.  Obviously,  the  basis 
multivectors  contain  the  basis  vectors  when  k  =  1. 

It  is  straightforward  to  generate  the  multiplication  table  of  the  basis  multivectors  of  the 
Clifford  algebra,  because  ej  ej  =  ±ex  for  some  compound  index  K .  The  procedure  is  as 
follows.  Move  each  factor  ej  for  q  =  1, . . . ,  l  in  the  product  c,jx  Cj2  ...  Cjk  Cjx  ej2  . . .  Cjt  to 
the  left  by  swapping  it  with  the  immediate  left  neighbour  ejp  while  ip  >  jq ,  and  changing 
the  sign  of  the  product  according  to  the  anti-commutativity  law.  Stop  when  either  ip  <  jq 
or  ip  =  jq.  In  the  latter  case  set  ejp  ejq  =  1.  In  the  end  of  the  procedure  the  indices  of 
the  product  factors  become  sorted  which  produce  the  final  index  K  and  the  sign  of  the 
equivalent  product.  For  example, 

e23  ei2  =  e2  es  ei  e2  =  —  e2  ei  e 3  e2  =  ei  e2  e 3  e 2  =  —e\  C2  C2  cz  =  —c\  ez  =  —  ei3  .  (6) 

In  particular,  we  have  proved  that  a  general  multivector  A  is  expressed  by  the  formal  sum 

A  =  YJMei  (7) 

1 

where  Aj  are  real  numbers. 

The  multiplication  rule  is  extended  to  all  multivectors  by  linearity  using  the  distribu¬ 
tivity  law.  By  virtue  of  the  formal  sum  (7),  we  immediately  arrive  at  the  direct  sum 
decomposition  of  the  Clifford  algebra 

C  (Mn)  =  C0  (Kn)  ©  Ci  (Mn)  ©  ...  ©  Cn  (Mn)  (8) 
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where  Ck  (Mn)  is  spanned  by  those  ej  whose  index  I  is  composed  of  k  integers.  The 
number  k  is  called  the  grade  of  the  vector  subspace  C j.  (Rn),  and  the  elements  of  C y.  (Mn) 
are  called  fc-vectors.  The  dimension  of  Ck  (Rn)  is  equal  to  ‘n  choose  k\  namely 

d'mC‘(R")  =  Q"U^U'  ,9) 

In  particular, 


dirnC  (Mn)  =  dim  Ck  (Rn)  = 

k= 0  k= 0 

Hence  the  Clifford  algebra  generated  by  an  n-dimensional  vector  space  is  2n-dimensional. 

The  vector  subspace  Co  (Mn)  is  isomorphic  to  M,  which  is  also  a  central  subalgebra. 
By  definition,  a  subalgebra  is  central  if  all  its  elements  commute  with  the  elements  of  the 
algebra  [Lang  2002],  The  centre  of  an  algebra  is  the  set  of  elements  commuting  with  all 
elements,  which  is  obviously  a  (maximal)  central  subalgebra  [Lang  2002].  It  is  known  that 
the  centre  of  C  (Mn)  is  Co  (Rn)  for  n  even,  or  Co  (Mn)  ©  Cn  (Mn)  for  n  odd. 

Note  that  0-vectors  are  called  scalars,  n-vectors  are  pseudoscalars,  and  1-vectors  C\  (Mn) 
are  naturally  identified  with  vectors  Rn.  The  following  direct  sum  of  the  even- grade  sub¬ 
spaces 

C+  (Mn)  =  C0  (Mn)  0  C2  (Rn)  0  ...  ©  C2Ln/2j  (Mn)  (11) 

is  in  fact  a  subalgebra  because  it  is  closed  under  multiplication.  Indeed,  the  product  of 
two  even-grade  elements  has  an  even  grade.  The  dimension  of  C+  (Mn)  is  equal  to  2n~1. 
The  even  subalgebra  C+  (Rn)  plays  an  important  role  in  the  definition  of  spinors  [Hestenes 
&  Sobczyk  1987,  Hurley  &  Vandyck  2000]. 

It  is  worthwhile  noting  that,  since  scalars  are  embedded  in  a  Clifford  algebra,  the 
multiplication  by  scalars  coincides  with  the  ordinary  multiplication  of  algebra  elements. 
Therefore,  from  the  algebraic  point  of  view,  Clifford  algebras  are  indeed  rings  [Lang  2002], 

Let  us  extend  the  inner  and  outer  products  from  1-vectors  to  all  multivectors.  A 
general  multivector  A  has  the  direct  sum  decomposition 

n 

n  =  (12) 

k= 0 

where  the  grade-projection  operator  of  C  (Mn)  to  Ck  (Mn)  is  denoted  by  the  angular  bracket 

C  (Mn)  (A)k  €  Ck  (Mn)  •  (13) 

It  is  convenient  to  extend  the  grade-projection  operator  to  identical  zero  for  k  >  n  or 
k  <  0.  It  suffices  to  define  the  inner  and  outer  products  on  elements  of  given  grades  and 
extend  them  by  linearity  to  all  multivectors.  Having  said  that  we  define 

n  n 

A-B  =  (14) 

k= 0 Z=0 


UNCLASSIFIED 


3 


DSTO-TR  3021 


UNCLASSIFIED 


^s  =  EE(w»Ww  <15) 

k= 0  1=0 

It  is  easy  to  verify  that  these  definitions  are  consistent  with  the  original  definitions  for 
1-vectors.  Note  that  a  -  A  =  A  -  a  =  a.  A  A  =  A  Aa  =  a  A  =  Aa  for  any  scalar  a  6  Co  (Mn) . 
The  inner  product  combines  the  left  and  right  contraction  in  a  single  expression.  It  is 
worthwhile  emphasizing  that  neither  is  the  inner  product  symmetric  nor  is  the  outer 
product  antisymmetric.  It  is  important  to  stress  that  the  outer  product  is  associative, 
whereas  the  inner  product  is  not. 

The  following  identities  exploiting  the  outer  and  inner  products  can  be  proved  in  a 


straightforward  manner 

a  ■  (b  A  c)  =  (a  ■  b)  c  —  (a  ■  c)  b  V  a,  b,  c  G  Mn  ,  (16) 

(a  A  b)  •  (c  A  d)  =  (a  •  d )  (b  ■  c)  —  (a  •  c)  (b  ■  d)  V  a,  b,  c,  d  €  Mn  ,  (17) 

aACAb  =  -bACAa  VC  6  C  (IT),  (18) 

a\  A  ...  A  an  =  det  [a\  ...  an]  e\  ...n  Vaj  G  Mn  .  (19) 


In  the  last  identity  the  expression  [ai  ...  an\  denotes  a  square  matrix  of  order  n,  formed 
by  the  components  of  the  n  vectors  arranged  as  columns  (or  rows). 

These  identities  can  be  generalized,  for  example,  the  identity  (16)  is  a  particular  case 
of  the  following  identity 

k 

a  •  B  =  ^(-l)*-1  (a  •  bi)  Bi  Vo,  &i,  . . . ,  bk  €  Mn  (20) 

i— 1 

where  B  =  b\  A  ...  A  bk  and  Bi  is  obtained  from  B  by  omitting  the  factor  6*.  Obviously, 
by  virtue  of  anticommutativity,  we  have  B  =  (— 1)*_1  bi  A  Bt. 

Another  important  operation  on  a  multivector  A,  denoted  by  A,  is  called  the  reversion. 
To  define  it,  express  a  multivector  A  in  the  form  (12)  and  reverse  the  order  of  the  factors 
in  the  basis  multivectors.  This  procedure  gives  the  expression 

n 

~  x — "\  fc  (fc— 1) 

A  =  (A),  .  (21) 

k=0 

In  particular,  the  reversion  has  no  effect  on  scalar  or  vector  components. 

In  order  to  represent  a  /c-dimensional  subspace  of  Mn  and  associate  algebraic  operations 
of  projection  and  rejection,  it  is  convenient  to  define  a  /c-blade  B  such  that  the  subspace 
becomes  the  kernel  of  the  map 

Rn3a^aABeC  (Mn)  .  (22) 

By  definition,  a  A:-blade  B  is  the  product  of  k  orthonormal  vectors,  namely 

B  =  h  . . .  bk  .  (23) 
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We  always  assume  that  k  >  1.  It  is  clear  that  the  blade  B  has  a  geometric  inverse  B  1 
given  by  the  reversion  B,  which  can  be  also  expressed  as 

B-1  =  bk  . . .  h  .  (24) 

Indeed,  B  B =  B~l  B  =  1  due  to  the  orthonormality  of  the  blade  factors.  Now,  if  the 
vectors  {b\,  . . . ,  b *.}  span  a  subspace  Mfc  C  Mn,  then  we  immediately  calculate  the  projected 
(parallel)  and  rejected  (perpendicular)  components  of  a  vector  a  by  the  following  formulae 

ci||  =  (a  ■  B)  B -1 ,  a±  =  (a  A  B)  B~1 .  (25) 

It  is  straightforward  to  verify  that  a  =  ay  +  a±  and  oy  •  a±  =  0.  These  simple  formulae 
have  no  analogue  in  Vector  Algebra  for  k  >  1.  In  particular,  for  a  vector  a  and  a  blade  B , 
the  identities  aB  =  a-  B  +  a/\B  and  B  a  =  B  ■  a  +  B  Aa  take  place,  which  generalize  the 
identity  ab  =  a ■ b  +  a  Ab  directly  following  from  the  definitions  (4) . 

Since  the  dimensions  of  C ^  (Rn)  and  Cra_£  (Mn)  are  equal,  there  exists  a  natural  isomor¬ 
phism  called  the  dual  operator,  which  is  closely  related  to  the  Hodge  star  operator.  The 
dual  operator  is  defined  by  the  following  map  of  basis  /c-vectors  to  basis  (n  —  fc)-vectors 

*ehi2  —ik  =  (— 1)1  I  ejij2  •••  jn-k 


where  the  compound  index  J  =  {j\,  j 2,  . . . ,  jn-k}  is  uniquely  determined  by  the  condition 
J  =  {1,  2,  . . . ,  n}  \  I  with  I  =  {?'i,  12,  . . . ,  ik},  and  |oj  is  the  parity  of  the  permutation 


a  = 


12  ...  n 

hh  ■  ■  ■  ik  jl  32  •  •  •  jn-k 


(27) 


Recall  that  1  <  i\  <  i^  <  ■  ■  ■  <  in  <  n  and  1  <  j\  <  ,72  <  •  •  •  <  jn-k  <  n-  In  particular, 
we  immediately  obtain  that  =t=l  =  e\2...n  and  *e\2...n  =  1-  It  is  straightforward  to  verify 
that  the  inverse  dual  operator  on  basis  multivectors  has  the  expression 


-1, 


-*1  *2  •••  Ik 


=  (-l 


M 


’'Jl  J2  •••  Jn  —  k 


(28) 


where  |r|  is  the  parity  of  the  permutation 


r  = 


12  ...  n 

.31  J2  •  •  •  jn-k  i\i2  ■  ■  ■  ik 


(29) 


The  dual  operator  (and  its  inverse)  is  extended  from  the  basis  multivectors  to  the  entire 
Clifford  algebra  by  linearity. 

Note  that  if  A  is  a  /c-blade,  1  <  k  <  n,  then  B  =  *A  is  a  (n  —  fc)-blade  representing  the 
orthogonal  complement  of  the  fc-dimensional  vector  subspace  defined  by  A.  In  particular, 
the  product  AB  is  a  nonzero  pseudoscalar  or,  equivalently,  a  n-blade. 

Let  us  define  the  anti- wedge  product  V  in  terms  of  the  wedge  product  A  and  the  dual 
operator  by  the  formula 


A  V  B  =  *  1  (*A  A  *B )  . 


(30) 
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Since  the  product  A  increases  and  V  decreases  grades,  the  wedge  and  anti-wedge  products 
are  also  called  the  progressive  and  regressive  outer  products,  respectively.  It  is  important 
to  stress  that  the  regressive  outer  product  is  also  associative.  Indeed,  since 

{AM  B)\J  C  =  *“1  (*A  /\*B  A  *C)  =  AM  (BV  C)  ,  (31) 

the  associativity  of  the  regressive  outer  product  is  the  direct  consequence  of  the  associa¬ 
tivity  of  the  progressive  outer  product. 

If  we  keep  only  the  progressive  outer  product  A  in  a  Clifford  algebra  but  all  the  other 
products  are  forgotten,  the  resulting  algebraic  structure  is  called  a  Grassmann  algebra. 
The  Grassmann  algebra  is  associative,  and  can  be  constructed  with  no  explicit  use  of  a 
scalar  product  in  the  original  vector  space.  For  example,  the  outer  product  can  be  defined 
in  terms  of  an  alternating  tensor  product.  Note  that  the  dual  isomorphism  in  a  Clifford 
algebra  relies  on  the  existence  of  a  non-degenerate  scalar  product  (not  necessarily  positive 
definite).  In  Grassmann  algebra,  which  has  less  structure  because  the  scalar  product 
may  not  naturally  exist,  a  similar  duality  can  be  established  by  specifying  a  nonzero 
pseudoscalar,  a  volume  element,  which  particularly  defines  an  orientation  of  the  vector 
space.  In  essence,  this  leads  to  the  definition  of  the  Grassmann-Cayley  algebra  which  is 
equipped  with  the  progressive  and  regressive  outer  products. 


2.1  Clifford  algebra  of  3-dimensional  vector  space 

Consider  a  3-dimensional  vector  space  R3  and  generate  the  corresponding  Clifford  algebra 
C  (R3)  spanned  by  8  basis  multivectors.  The  multiplication  table  of  the  basis  multivectors 
is  presented  in  Table  Al,  and  the  progressive  outer  products  are  given  in  Table  A2. 

The  dual  operator  is  calculated  on  the  basis  multivectors  in  a  straightforward  manner, 
namely 

*1  =  ei23  ,  *ei23  =  1 ,  (32) 

*ei  =  e23  ,  *e2  =  —  ei3  ,  *e3  =  ei2  ,  (33) 

*ei2  =  e3  ,  *ei3  =  — e2  ,  *e23  =  ■  (34) 

In  particular,  we  obtain  the  expression:  =  *. 

For  any  two  1-vectors  a,  b  we  get  the  identity  *  (a  A  b)  =  a  x  b  where  the  symbol  x 
denotes  the  cross  product  in  a  3-dimensional  Euclidean  space.  Note  that  the  cross  product 
is  not  associative,  whereas  the  outer  product  is.  The  dual  operator  allows  one  to  calculate 
the  regressive  outer  products  of  the  basis  multivectors.  The  result  is  shown  in  Table  A3. 

It  is  seen  that  the  scalar  unit  1  and  pseudoscalar  unit  ei23  commute  with  all  the 
other  basis  multivectors  with  respect  to  the  geometric  product,  and  therefore  the  subspace 
spanned  by  {l,ei23}  forms  the  central  subalgebra  of  the  Clifford  algebra  C  (R3),  which 
is  isomorphic  to  the  set  of  complex  numbers  C.  Actually,  this  central  subalgebra  is  the 
centre  of  the  algebra.  The  pseudoscalar  unit  ei23  plays  the  role  of  the  imaginary  unit. 
In  general,  any  element  a  6  C  (R3)  squaring  to  —1  (for  example  a  =  ei2)  generates  a 
2-dimensional  central  subalgebra  spanned  by  {1,  a},  which  is  isomorphic  to  C.  However, 
such  a  subalgebra  is  not  central. 
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The  even  subalgebra  C+  (R3)  is  isomorphic  to  the  algebra  of  quaternions  H  spanned 
by  the  basis  {1,  i,  j,  k}  satisfying  the  identities 

i2=j2  =  k2  =  — 1,  i  j  =  k  =  — j  i ,  j  k  =  i  =  — k  j  ,  ki=j  =  -ik.  (35) 

The  isomorphism  is  established  by  the  correspondence 

i  =  —623  =  -*e1:  j  =  ei3  =  -  *  e2  ,  k  =  -ei2  =  -  *  e3  .  (36) 

The  inverse  of  a  nonzero  quaternion  q  E  C+  (R3)  is  given  by  the  formula 

q -1  =  4  (37) 

qq 

where  the  product  q  q  is  obviously  a  positive  scalar.  The  quaternion  q  induces  a  rotation 
of  the  vector  space  R3  by  the  map  a  H >  pq(a )  =  qaq _1.  Indeed,  it  is  straightforward  to 
demonstrate  that  pq(a)  E  R3  and  pq(a)  ■  pq(b)  =  a  ■  b  for  all  a,  b  E  R3.  In  addition,  this 
map  preserves  the  space  orientation:  pq(c)  =  pq(a)  x  pq(b)  whenever  c  =  a  x  b. 

Note  that  the  quaternion  q  is  determined  by  the  rotation  pq  up  to  a  nonzero  scale. 
In  particular,  the  multiplicative  group  of  normalized  quaternions,  which  is  a  unit  sphere 
{qq  =  1}  in  the  four-dimensional  space  of  quaternion  components,  is  a  double  covering  of 
the  group  of  rotations  of  a  three-dimensional  Euclidean  space  (pq  =  p~q).  The  map  q  pq 
of  the  group  of  normalized  quaternions  to  the  group  of  rotations  is  a  group  homomorphism 
since  pqiq2  =  pqi  pq2  [Lang  2002], 

There  are  several  other  subalgebras  such  as  that  spanned  by  {1,  ei,  e2,  ei2},  which  is 
isomorphic  to  C  (R2).  In  general,  any  pair  of  orthonormal  vectors  generates  a  subalgebra 
of  C  (R3) ,  which  is  isomorphic  to  C  (R2) . 

The  less  trivial  example  is  the  commutative  subalgebra  spanned  by  {1,  ei,  e23,  ei23}, 
which  corresponds  to  tessarines  also  called  bicomplex  numbers.  In  general,  any  vector  a 
squaring  to  1  generates  a  tessarine  subalgebra  spanned  by  {1,  a,  a  ei23,  ei23}. 


2.2  Clifford  algebra  of  4-dimensional  vector  space 

Consider  a  4-dimensional  vector  space  R4  and  generate  the  corresponding  Clifford  algebra 
C  (R4)  spanned  by  16  basis  multivectors.  The  multiplication  table  of  the  basis  multivectors 
is  presented  in  Table  A4,  and  the  progressive  outer  products  are  given  in  Table  A5.  The 
dual  operator  is  calculated  on  the  basis  multivectors  as  follows 


*1  =  ei234  , 

*ei234  =  1 , 

(38) 

*ei  =  e234 , 

*e2  =  — ei34  , 

*e3  =  ei24  , 

*e4  =  — ei23  , 

(39) 

*ei23  =  e4 , 

*ei24  =  —  e3  , 

*ei34  =  e2  , 

*e234  =  — ei , 

(40) 

*ei2  =  e34 , 

*ei3  =  -e24 , 

*ei4  =  e23  , 

(41) 

*e23  =  ei4 , 

*e24  =  — ei3  , 

=t=e34  =  ei2  • 

(42) 
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Note  that  *_1  7^  *.  The  1-vectors  are  dual  to  3- vectors,  whereas  the  2- vectors  are  self-dual. 
There  are  many  subalgebras  of  C  (R4),  such  as  C  (R3)  generated  by  any  three  orthonormal 
vectors  of  R4.  In  particular,  Table  A1  is  a  subtable  of  Table  A4. 

This  algebra  naturally  arises  in  Projective  Geometry  that  models  the  3-dimensional  Eu¬ 
clidean  space  completed  with  a  plane  at  infinity  using  four  homogeneous  coordinates.  Note 
that  in  Physics  a  4-dinrensional  vector  space  equipped  with  an  indefinite  scalar  product 
with  signature  (+,  — ,  — ,  — )  models  Minkowski  spacetime,  so  the  multiplication  tables  and 
the  dual  operator  are  computed  differently.  Alternatively,  the  Minkowski  spacetime  can 
be  identified  with  the  subspace  Cq  (R3)  ©  Ci  (R3)  whose  elements  are  called  paravectors. 


3  Projective  geometry 

Consider  an  n-dimensional  projective  space  Pra  over  real  numbers  R,  which  is  identified 
with  the  set  of  all  straight  lines  through  the  origin  in  Rn+1  [Casse  2006].  It  is  convenient  to 
introduce  homogeneous  coordinates  (aq,  X2,  ■  ■  ■ ,  xn+  \ )  specifying  a  particular  line  in  Rn+1, 
a  point  of  Pn,  assuming  that  at  least  one  of  the  Xi  is  nonzero.  A  point  of  Pn  is  defined  by 
(xi,X2,  •  •  • ,  xn+\)  up  to  a  nonzero  scale.  Even  if  we  normalize  the  homogeneous  coordinates 
by  the  condition 

xj  +  x\  +  •  •  •  +  xn+ 1  =  1  ?  (43) 

the  corresponding  line  will  still  be  defined  by  (xi,X2,  ■  •  • ,  xn+i)  up  to  a  factor  ±1.  So,  from 
the  topological  point  of  view,  a  projective  space  P”  is  obtained  from  an  n-dimensional 
sphere  defined  by  (43)  by  identifying  (glueing)  diametrically  opposite  points.  For  a  point 
X  G  P"  we  denote  homogeneous  coordinates  by  X  =  (aq,X2, . . . ,  a;n+i)  G  Rn+1  keeping  in 
mind  the  unimportance  of  an  overall  scaling. 

It  is  worthwhile  noting  that  the  projective  space  Pn  includes  Rn  through  the  following 
map  called  homogenization 

Rn  9  (xi,x2,  ■  .,jXn)  |-t  {xi,X2,  ...,xn,l)  G  Rn+1 .  (44) 

The  inverse  projection  called  dehomogenization  is  defined  by  the  map 

Rn+1  9  (xi,x2,  •  •  .,xn+i)  i-A  -^2-,  •  •  • ,  G  Rn  (45) 

V^n+l  ^n+1  ^n+1 J 

which  makes  sense  only  if  xn+ 1  7^  0.  By  definition,  a  point  of  P"  lies  at  infinity  if  xn+i  =  0. 

Construct  the  Clifford  algebra  C  (Rn+1)  using  the  set  of  homogeneous  coordinates  as 
the  underlying  vector  space.  Note  that  two  points  X\,  X2  of  P"  are  equal  if  and  only 
if  X[  A  X-2  =  0.  Let  X\,  X2  be  two  distinct  points  of  Pn  represented  by  homogeneous 
coordinates  X\ .  X2.  In  particular,  the  1-vectors  X\ ,  X2  are  not  collinear.  Define  the 
nonzero  2-vector  L  =  X\  A  X2  which  represents  a  two-dimensional  plane  in  Rn+1  passing 
though  X\ ,  X2  and  the  origin.  The  dehomogenization  (45)  maps  this  plane  to  a  line  L 
in  Rn  passing  through  X\  and  A2.  If  both  X\  and  X2  are  at  infinity,  the  line  L  will  be 
also  at  infinity,  therefore  L  is  a  line  in  P".  This  line  as  a  set  is  defined  by  the  kernel  of 
the  map 

Rn+1  Al  AieC3  (Rn+1)  .  (46) 
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Clearly,  X\  A  L  =  X2  A  L  =  0  or,  equivalently,  X\,  X2  €  L  by  construction.  Note  that  the 
overall  scales  of  the  homogeneous  coordinates  do  not  change  the  line  as  a  set. 

Likewise,  three  points  Xi,  X2,  X3  of  a  projective  space  Pn  whose  vectors  of  homo¬ 
geneous  coordinates  are  not  coplanar  with  the  origin  of  Mn+1  define  a  two-dimensional 
plane  P  in  Pn  by  the  nonzero  3-vector  P  =  X\  A  X2  A  X3  (we  assume  n  >  2).  This  plane 
as  a  set  is  defined  by  the  kernel  of  the  map 

Mn+1  3X4XA  PeC4  (Mn+1)  .  (47) 

Clearly,  X\  A  P  =  X2  A  P  =  X3  A  P  =  0,  thus  reinstating  the  fact  that  X \ ,  X2 ,  A3  €  P  by 
construction.  A  plane  P  can  be  also  obtained  from  a  line  L  represented  by  a  2-vector  L 
and  a  point  X  represented  by  a  1- vector  X  using  the  product  P  =  L  A  X. 

In  other  words,  the  progressive  outer  product  provides  an  efficient  way  to  build  objects 
of  higher  dimension  from  objects  of  lower  dimension.  This  operation  is  called  a  ‘join’.  The 
complementary  operation  called  a  ‘meet’  is  handled  by  the  regressive  outer  product,  and 
corresponds  to  the  intersection  of  objects.  It  is  worthwhile  emphasizing  that,  though  the 
algebraic  operations  of  Clifford  algebra  are  always  carried  out,  the  result  can  be  zero.  This 
situation  indicates  that  the  resulting  object  is  invalid.  For  example,  this  happens  when 
two  identical  objects  are  multiplied. 


3.1  Incidence  relations  in  projective  plane  P2 

Consider  two  distinct  lines  Li  and  L2  in  P2  represented  by  two  2- vectors  in  C  (M3) .  Com¬ 
pute  the  regressive  product  X  =  L\  V  L2  which  is  obviously  a  1-vector  and  therefore 
represents  a  point  X  of  P2.  It  is  straightforward  to  prove  that  X  is  the  intersection  point 
of  the  lines  Li  and  L2  possibly  lying  at  infinity.  Also  note  that  a  point  X  belongs  to  a 
line  L  if  either  X  A  L  =  0  or  X  V  L  =  0. 

Two  points  X\ ,  X2  of  the  projective  plane  P2  are  equal  if  either  X\  A  X2  =  0  or 
X\  \J X2  =  0.  Likewise,  two  lines  L2  are  equal  if  either  Li  AL2  =  0  or  L1VL2  =  0.  There 
is  a  natural  duality  between  points  and  lines  established  by  the  mutual  correspondences 
L  =  *X  and  X  =  *L. 

Let  A\ .  B 1,  Ci  and  A2,  B2,  C2  be  two  triplets  of  collinear  points  in  the  projective 
plane  P2.  Algebraically,  this  means  that  A\  A  B\  A  C\  =  0  and  A2  A  B2  A  C2  =  0.  Then 
the  intersection  points  A,  B,  C  constructed  by  the  rules 

A=(b  1  a  c2)  V  (b2  a  Cl)  ,  (48) 

B  =  (Cy  A  i2)  V  (C2  A  ii)  ,  (49) 

C  =  ^Ai  A  -B2)  V  ^A_2  A  Bi'j  ,  (50) 

must  be  also  collinear,  that  isj4ABAC  =  0.  This  is  the  statement  of  the  theorem  of 
Pappus.  A  dual  version  of  this  theorem  exists  in  which  collinear  points  are  replaced  with 
concurrent  lines. 
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Now  assume  that  the  points  A\,  B i,  C\  and  A 2,  B2,  C2  are  not  collinear  but  form  two 
non-degenerate  triangles.  Construct  the  intersection  points  A,  B,  C  of  the  corresponding 
triangle  sides  by  the  rules 

A=(b1  a  Ci)  V  (b2  a  c2)  ,  (51) 

B  =  (dy  A  ii)  V  (62  A  i2)  ,  (52) 

C  =  (ii  A  Pi)  V  (i2  A  P2)  ,  (53) 

and  compute  the  pseudoscalar  p  =  A/\B  AC  which  vanishes  if  and  only  if  the  intersection 

points  are  collinear.  Draw  three  lines  La ,  Lb,  Lc  through  the  corresponding  vertices  of  the 
triangles,  namely 

La  =  A\  A  A2  ,  Lb  =  B\  A  P2  ,  Lc  =  C\  A  C*2  ,  (54) 

and  compute  the  scalar  q  =  La\l  Lb\!  Lc  which  vanishes  if  and  only  if  the  lines  are 
concurrent.  The  theorem  of  Desargues  states  that  p  and  q  vanish  simultaneously.  This 
statement  follows  from  the  explicit  formula 

*p  =  —I]  I2  q  ,  I\  =  A\  A  B\  A  C\  ,  I2  =  A2  A  -B2  A  C2  ,  (55) 

in  combination  with  the  fact  that  the  pseudoscalars  I±,  /2  do  not  vanish  (recall  that  the 
triangles  are  not  degenerate). 


3.2  Incidence  relations  in  projective  space  P3 


Now  consider  a  line  L  and  a  plane  P  in  P3  represented  respectively  by  an  appropriate 
2- vector  L  and  a  3-vector  P  in  C  (M4) .  Then  the  product  X  =  L  V  P,  which  is  obviously  a 
1-vector,  represents  the  intersection  point  X.  The  intersection  of  two  planes,  Pi  and  P2, 
is  given  by  the  product  L  =  Pi  V  P2  which  is  a  2-vector  in  C  (R4)  and  therefore  represents 
a  line  L.  The  point  of  intersection  of  three  planes  is  represented  by  the  triple  product 
X  =  P\  V  P2  V  P3  which  is  obviously  a  1-vector.  These  operations  work  for  general  (non¬ 
degenerate)  configurations. 

Using  either  the  regressive  or  progressive  product,  multiply  a  point  and  a  plane  to  get 
a  signed  minimum  distance  between  them.  Likewise,  multiply  two  lines  Li  and  L2  to  get  a 
special  signed  crossing  value.  Depending  on  the  product,  the  result  will  be  either  a  scalar 
or  pseudoscalar.  If  the  crossing  value  vanishes,  say  L\  A  P2  =  0,  the  lines  are  coplanar  and 
hence  intersect.  The  homogeneous  coordinates  of  the  intersection  point  X  belong  to  the 
kernel  of  the  map 

143I^(1aIi,1aL2)  €  C3  (R4)  x  C3  (R4)  .  (56) 

In  other  words,  we  have  to  find  a  nonzero  solution  to  the  system  of  equations 


(  X  A  Li  =  0 
\  X  A  L2  =  0 


(57) 
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which  is  equivalent  to  an  overdetermined  system  of  8  homogeneous  linear  equations  for  4 
components  of  X.  The  solution  can  be  efficiently  obtained  by  the  Singular  Value  De¬ 
composition  (SVD)  algorithm  of  Linear  Algebra  [Golub  &  Van  Loan  1996].  The  matrix 
equation  is  derived  in  a  straightforward  manner  by  observing  that  the  equation  X  A  L  =  0 
is  equivalent  to 


1 

to 

CO 

1 

CO 

to 

0 

_ 1 

’  Xl 

'  0  ' 

hi  ~hi  0  l\2 

X2 

0 

co'4 

O 

CO 

X3 

0 

|_  0  ^34  ~hi  I23 

X4 

0 

where 

X  =  x\  e\  +  x'2  e2  +  X3  e3  +  x4  e 4  , 

L  =  l\2  ei2  +  ^13  ei3  +  Z14  ei4  +  I23  e23  +  hi  ^24  +  ^34  634  . 


(58) 


(59) 

(60) 


Note  that  the  six  coefficients  in  the  line  representation  (60)  are  the  Pliicker  coordinates. 

The  Pliicker  coordinates  providing  an  economic  parametrization  of  all  lines  in  P3  are 
not  arbitrary  but  satisfy  a  constraint  imposed  by  the  Klein  quadric 


h-2  hi  ~  I13  hi  +  ^14  I23  —  0  . 

A  line  L  =  X  A  Y  constructed  by  two  points  has  the  expression 

L  =  (xi  y2  -  ®2?/i)  ei2  +  (xu/3  -  x3yi)  ei3  +  (xij/4  -  x4yi)  ei4 
+  (X2  V3  -  x3  2/2)  e23  +  (x2  Vi  -  X4  V2)  e24  +  (x3  y4  -  xA  y3)  e34 


(61) 


(62) 


which  no  longer  contains  any  information  about  the  1-vectors  X,  Y  used  to  create  it.  This 
is  contrary  to  a  parametric  representation  of  a  line  L  passing  through  points  X,  Y.  It  is 
straightforward  to  check  that  the  constraint  (61)  is  imposed  for  the  line  (62). 

Two  points  Xi,  V2  of  the  projective  space  P3  are  equal  if  X\  A  V2  =  0.  Likewise,  two 
planes  Pi,  P2  are  equal  if  P\  V  P2  =  0.  There  is  a  natural  duality  between  points  and 
planes  established  by  the  correspondences  P  =  *X  and  X  =  *P. 


The  incidence  relations  with  lines  are  more  involved.  Two  lines  Pi,  P2  are  equal  if  the 
kernel  of  the  map  (56)  is  two-dimensional.  This  question  can  be  efficiently  answered  using 
the  SVD  algorithm,  or  by  the  direct  calculation  of  the  rank  r  of  the  matrix 


M( Pi,P2) 
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7(2) 

‘12 

0 

7(2) 
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'34 

7(2) 

‘24 

7(2) 
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(63) 


where  are  the  Pliicker  coordinates  of  the  line  Lm.  The  lines  are  not  coplanar  if  r  =  4, 
the  lines  have  a  unique  point  of  intersection  if  r  =  3,  and  the  lines  coincide  if  r  =  2.  The 
case  r  <  2  never  happens  because  a  valid  line  as  a  set  is  a  one-dimensional  object. 
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4  Application  to  3D  reconstruction 

In  this  section  we  provide  a  theoretical  background  for  3D  reconstruction  using  Geometric 
Algebra.  We  derive  several  important  concepts  and  relations  from  [Hartley  &  Zisserman 
2003]  in  a  simple  and  elegant  form. 

The  three-dimensional  space  will  be  modelled  by  the  projective  space  P3.  To  apply 
efficient  algebraic  methods  when  building  geometric  objects,  we  generate  the  Clifford  alge¬ 
bra  C  (M4) .  We  will  slightly  abuse  notation  by  using  the  same  symbol  for  an  object  in  P3 
and  its  representation  in  C  (M4).  Recall  that  the  object  representation  in  C  (M4)  is  defined 
up  to  a  nonzero  scale. 


Figure  1:  Pinhole  camera 


A  pinhole  camera  is  completely  described  by  its  optical  centre  C  and  projection  plane  P 
both  considered  to  be  objects  of  the  projective  space  P3.  It  is  always  guaranteed  that  C 
does  not  belong  to  P,  which  is  equivalent  to  C  A  P  /  0.  The  projective  space  P3  has 
no  natural  metric  structure.  However,  image  points  are  measured  in  Cartesian  coordi¬ 
nates  of  the  projection  plane  which  does  have  a  metric.  So,  it  is  important  to  establish  a 
parametrization  of  image  points  in  P3  in  terms  of  the  Cartesian  coordinates  of  the  projec¬ 
tion  plane.  This  can  be  done,  for  example,  by  identifying  some  point  O,  say  the  principal 
point,  with  the  origin  of  the  coordinate  system  and  choosing  two  points  A,  B  on  the 
axes  of  the  orthogonal  coordinates.  The  four  points  C,  O ,  A,  B  completely  determine 
the  pinhole  camera,  in  particular  P  =  O  A  A  A  B.  A  point  X  in  the  projection  plane  P 
has  the  representation  X  =  a  A  +  b  B  +  cO  where  the  scalars  a,  b ,  c  are  defined  up  to  a 
common  scale.  These  scalars  are  directly  related  to  the  Cartesian  coordinates.  Indeed, 
since  all  the  points  X,  O,  A ,  B  are  finite  (not  lying  at  infinity),  we  can  dehomogenize  them 
and  set  a  +  b  +  c  =  1  without  loss  of  generality.  As  a  result,  we  get  the  representation 
X  =  O  +  a  {A  —  O)  +  b  (B  —  O)  where  a,  b  are  the  Cartesian  coordinates.  The  direction 
vectors  A  =  A  —  O  and  B  =  B  —  O  belong  to  the  plane  at  infinity  of  P3 ,  and  are  orthogonal 
by  construction.  The  Cartesian  coordinates  are  now  readily  computed  by  the  formulae 


a  = 


AX 
A- A’ 


b 


B  ■  X 
B  ■  B  ' 


(64) 
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where  X  =  X  —  O  is  the  direction  vector.  The  scalar  product  is  induced  by  the  projection 
plane  P. 

It  is  important  to  emphasize  that  it  makes  sense  to  talk  about  orthogonality  of  par¬ 
ticular  vectors  constructed  from  the  projection  plane  P.  In  fact,  most  problems  of  3D 
reconstruction  are  formulated  as  a  mathematical  problem  of  fitting  a  projective  (or  affine) 
variety  [Cox,  Little  &  O’Shea  2007]  to  measurement  points  by  minimizing  a  cost  function 
associated  with  the  metric  structure  of  projection  planes. 

Let  X  E  IP3  be  an  arbitrary  point  different  from  the  optical  centre  C.  The  following 
map  projects  the  point  X  to  the  projection  plane  P 

I^PV(CAI).  (65) 

This  map  generates  the  camera  matrix  when  the  projected  point  is  represented  in  the 
Cartesian  coordinates  of  the  projection  plane  P  followed  by  the  homogenization  procedure. 
The  size  of  the  camera  matrix  is  3-by-4  as  it  maps  four  homogeneous  coordinates  of  the 
point  X  to  the  three  homogeneous  coordinates  a,  b ,  c  of  the  projection  plane. 

Now  let  L  be  an  arbitrary  line  in  P3  not  containing  the  optical  centre  C.  The  projection 
of  this  line  to  the  plane  P  is  given  by  the  similar  map 

L  H-  P  V  (C  A  I)  .  (66) 

The  result  is  the  unique  line  of  intersection  of  two  planes. 


Figure  2:  Epipolar  geometry 

In  multiview  geometry  we  are  provided  with  a  sequence  of  optical  centres  and  projection 
planes  which  we  label  by  some  index.  Given  two  camera  configurations  ( CQ,Pa )  and 
(Cp,  Pp)  we  construct  the  base  line  Ca/\Cp  and  intersect  it  with  the  projection  planes  (see 
Figure  2).  The  intersection  points  are  respectively  the  epipoles  Oap  and  Opa.  The  epipolar 
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lines  Eap  and  Epa  are  formed  by  the  intersection  of  a  plane  containing  the  baseline  with 
the  projection  planes.  Obviously,  all  the  epipolar  lines  in  a  projection  plane  are  concurrent 
as  they  contain  the  epipole,  and  therefore  have  the  geometry  of  the  projective  line  P1. 

There  exists  the  fundamental  mapping  <&ap  between  the  epipolar  lines,  induced  by  the 
pencil  of  planes  containing  the  baseline.  Algebraically,  the  fundamental  mapping  is  given 
by  the  formula 


<F 


'a/3  '■  IP1  3  Eaf)  i— >•  Epa  —  (Ca  A  Eap)  V  Pp  G 


(67) 


It  is  straightforward  to  prove  that  the  mapping  &ap  is  well  defined  (recall  that  Cp  £  Pp) 
and  is  invertible.  Indeed,  by  virtue  of  symmetry,  we  have  Qpa  = 

It  is  clear  that,  if  Xa  G  Pa  and  Xp  G  Pp  are  two  images  of  some  point  X,  then 
Xa  G  Eap  and  Xp  G  Epa  as  shown  in  Figure  2.  Algebraically,  this  is  equivalent  to  the 


constraint 


$a/3  {Xa  a  Oap)  A  Xp  =  0 


(68) 


where  the  epipole  has  the  expression  Oap  =  Pa  V  {Ca/\Cp).  This  equation  produces 
the  fundamental  matrix  [Hartley  &  Zisserman  2003]  when  Xa  and  Xp  are  respectively 
represented  in  coordinates  of  the  planes  Pa  and  Pp. 

Given  two  image  points  Xa,  Xp  satisfying  the  constraint  (68),  the  point  X  is  recovered 
by  the  intersection  of  the  optical  rays  CQ  A  XQ  and  Cp  A  Xp.  This  operation  called 
triangulation  is  well  defined  if  the  optical  rays  are  coplanar.  The  later  is  equivalent  to  (68) 
but  a  simpler  form  which  does  not  involve  the  projection  plane  Pp  is  as  follows 


Xa  A  Ca  A  Cp  A  Xp  =  0  . 

The  triangulation  problem  reduces  to  the  system  of  equations 

X  A  CQ  A  Xa  =  0 
X  A  Cp  A  Xp  =  0 


(69) 


(70) 


which  is  of  the  form  (57).  A  nonzero  solution  (up  to  a  scaling  factor)  exists  provided  that 
the  constraint  (69)  is  imposed.  The  conditions  for  Xa  G  Pa  and  Xp  G  Pp  are  equivalent 
to  the  following 


Xa  A  Pa  =  0  ,  Xp  A  Pp  =  0  . 


(71) 


Here  the  progressive  product  in  the  inclusion  tests  can  be  replaced  with  the  regressive 
product,  since  points  and  planes  are  in  the  duality  relationship. 

It  is  worthwhile  emphasizing  that  the  projection  planes  do  not  appear  explicitly  in 
the  triangulation  problem.  This  is  because  we  consider  the  image  points  belonging  to  the 
projective  space  P3. 

The  triangulation  procedure  solves  the  following  problem  of  3D  reconstruction  in  a 
closed  form:  “Given  two  camera  configurations  and  arrays  of  pairs  of  matching  image 
points  satisfying  the  constraint  (69)  exactly,  find  the  corresponding  spatial  points.” 

However,  in  reality,  the  3D  reconstruction  procedure  is  complicated  by  the  presence  of 
unavoidable  noise  in  the  measurements  of  image  points.  In  particular,  we  cannot  assume 
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any  longer  that  the  constraint  (69)  is  exactly  satisfied.  Another  problem  is  associated  with 
a  priori  unknown  camera  configurations.  Finally,  the  matching  points  when  automatically 
produced  may  not  be  in  correspondence  at  all.  Such  points  are  called  outliers  as  opposed 
to  inliers,  and  must  be  eliminated.  The  fundamental  matrix  plays  the  crucial  role  in  the 
elimination  procedure  [Hartley  &;  Zisserman  2003]. 

In  the  two- view  geometry  the  reconstruction  of  lines  is  less  restrictive.  Let  La  and  Lg 
be  two  image  lines  in  the  planes  Pa  and  Pg ,  not  containing  the  epipoles  Oag  and  Oga, 
respectively.  In  other  words,  neither  of  them  is  an  epipolar  line.  Then  the  reconstructed 
line  L  is  given  by  the  formula 

L  =  ( Ca  A  La )  V  (Cg  A  Lg)  .  (72) 

The  result  is  the  unique  line  of  intersection  of  two  different  planes.  However,  this  advantage 
turns  out  to  be  disadvantage  as  there  is  no  way  to  eliminate  outlying  lines.  The  situation 
becomes  different  in  the  three- view  geometry  [Hartley  &:  Zisserman  2003]. 

If  both  projection  lines  are  epipolar  lines,  then  the  planes  are  either  equal  or  intersect 
at  the  base  line.  The  base  line  is  obviously  projected  to  the  epipoles  of  the  projection 
planes,  hence  this  situation  cannot  generate  image  lines.  If  only  one  projection  line  is 
epipolar,  say  La,  but  the  other  line  Lg  is  not,  then  both  planes  Ca  A  La  and  Cg  A  Lg 
contain  the  optical  centre  Cg.  Therefore,  the  reconstructed  line  L  contains  Cg  and  must 
be  projected  to  the  epipole  OgQ.  This  can  be  reformulated  in  the  following  way.  If  a  line  L 
contains  the  optical  centre  of  one  view,  then  in  the  other  view  it  must  be  projected  to 
either  an  epipolar  line  or  the  epipole. 

As  in  the  triangulation  problem,  the  projection  planes  do  not  appear  explicitly  in  the 
solution.  The  conditions  for  La  C  Pa  and  Lg  C  Pg  are  equivalent  to  the  following 

La  V  Pa  =  0  ,  Lg  V  Pg  =  0  .  (73) 

Here  the  regressive  product  in  the  inclusion  tests  cannot  be  replaced  with  the  progressive 
product. 


5  Description  of  the  developed  MATLAB  code 

In  this  section  we  briefly  describe  the  MATLAB  code  developed.  Basic  algebraic  operations 
of  the  Clifford  algebra,  currently  supporting  dimensions  up  to  n  =  4,  are  implemented  in 
the  cliff ord_algebra.m  file  whose  content  is  provided  in  Appendix  B.  This  reusable 
MATLAB  function  has  to  be  placed  on  the  MATLAB  path  to  be  called  from  different 
directories. 

For  the  sake  of  code  performance,  we  have  deliberately  avoided  the  construction  of 
MATLAB  classes.  Instead,  a  collection  of  multivectors  is  represented  by  an  ordinary 
matrix  containing  2n  columns.  The  fully  vectorized  code  is  capable  of  operating  on  large 
arrays  of  multivectors  in  an  efficient  way.  Factors  in  binary  products  can  be  either  arrays  of 
multivectors  of  the  same  number  of  rows  or  an  array  and  a  single  multivector  represented 
by  a  row  matrix. 
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The  unit  tests  of  the  Clifford  algebra  implementation  are  provided  in  Appendix  C. 
The  basic  properties  of  algebraic  operations  and  known  identities  derived  in  Clifford  alge¬ 
bra  are  thoroughly  tested  using  randomly  generated  arrays  of  multivectors.  For  example, 
we  have  validated  the  associativity  of  geometric,  progressive  and  regressive  products,  the 
basic  properties  of  scalar,  inner  and  cross  products,  and  the  properties  of  subalgebras 
including  even  subalgebra,  complex  and  bicomplex  numbers,  and  quaternions.  The  intro¬ 
duced  operations,  such  as  the  dual,  reverse  and  blade  operations,  are  also  tested.  The 
last  test  displays  the  performance  of  the  vectorized  geometric  product  against  its  serial 
counterpart. 

The  application  of  Clifford  algebra  to  Projective  Geometry  is  given  in  Appendix  D. 
The  unit  tests  include  the  validation  of  basic  incidence  relations  and  celebrated  theorems 
of  Projective  Geometry  such  as  those  of  Pappus  and  Desargues.  Also,  epipolar  geometry, 
fundamental  map,  and  the  point  and  line  reconstruction  procedures  are  tested. 


6  Discussion 

We  have  presented  a  fully  vectorized  MATLAB  code  for  the  basic  operations  of  Geometric 
Algebra  whose  real  power  manifests  itself  when  particularly  combined  with  Projective  Ge¬ 
ometry.  The  developed  code  is  thoroughly  tested  and  can  be  used  in  various  applications 
which  involve  complicated  geometric  constructs.  The  vectorized  code  can  be  applied  to 
real-time  simulations  when  there  is  a  requirement  to  operate  on  large  arrays  of  multivec¬ 
tors.  For  example,  quaternions  naturally  embedded  in  the  Clifford  algebra  C  (M3)  can  be 
utilized  to  implement  fast  rotation  of  large  arrays  of  3D  objects. 

Having  mainly  in  mind  the  application  of  the  Clifford  algebra  to  the  reconstruction  of  a 
three-dimensional  structure  from  a  sequence  of  images,  we  confined  ourselves  to  a  positive 
definite  scalar  product  when  defining  the  algebra  multiplication  rules.  The  extension  of 
the  code  to  an  indefinite  scalar  product  which  arises  in  the  geometry  of  spacetime  is 
straight  forward . 
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Appendix  A  Multiplication  tables 

The  multiplication  tables  of  two  Clifford  algebras,  C  (M3)  and  C  (M4),  are  given  in  this  ap¬ 
pendix.  The  basis  multivectors  in  the  leftmost  column  and  in  the  top  row  are  respectively 
the  left  and  right  factors  (multiplicand  and  multiplier)  of  the  binary  products. 


Table  Al:  Geometric  products  of  the  basis  multivectors  of  C  (M3) 


1 

ei 

e2 

e3 

ei2 

ei3 

e23 

ei23 

1 

1 

ei 

e2 

e3 

ei2 

ei3 

e23 

ei23 

ei 

ei 

1 

ei2 

ei3 

62 

63 

ei23 

e23 

e2 

e2 

— ei2 

1 

e23 

-ei 

— ei23 

63 

— ei3 

e3 

e3 

— ei3 

— e23 

1 

ei23 

-ei 

—  62 

ei2 

ei2 

ei2 

—e-2 

ei 

ei23 

-1 

— e23 

ei3 

-e3 

ei3 

ei3 

-e3 

— ei23 

ei 

e23 

-1 

— ei2 

e2 

e23 

e23 

ei23 

-e3 

e2 

— ei3 

ei2 

-1 

-ei 

ei23 

ei23 

e23 

— ei3 

ei2 

—63 

62 

-ei 

-1 

Table  A2:  Progressive  outer  products  of  the  basis  multivectors  of  C  (M3) 


A 

1 

ei 

e2 

e3 

ei2 

ei3 

e23 

ei23 

1 

1 

ei 

e2 

e3 

ei2 

ei3 

e23 

ei23 

ei 

ei 

0 

ei2 

ei3 

0 

0 

ei23 

0 

e2 

e2 

— ei2 

0 

e23 

0 

— ei23 

0 

0 

e3 

e3 

— ei3 

— e23 

0 

ei23 

0 

0 

0 

ei2 

ei2 

0 

0 

ei23 

0 

0 

0 

0 

ei3 

ei3 

0 

— ei23 

0 

0 

0 

0 

0 

e23 

e23 

ei23 

0 

0 

0 

0 

0 

0 

ei23 

ei23 

0 

0 

0 

0 

0 

0 

0 

Table  A3:  Regressive  outer  products  of  the  basis  multivectors  of  C  (M3) 


V 

1 

ei 

e2 

e3 

ei2 

ei3 

e23 

ei23 

1 

0 

0 

0 

0 

0 

0 

0 

1 

ei 

0 

0 

0 

0 

0 

0 

1 

ei 

e2 

0 

0 

0 

0 

0 

-1 

0 

e2 

0 

0 

0 

0 

1 

0 

0 

e3 

ei2 

0 

0 

0 

1 

0 

ei 

e2 

ei2 

ei3 

0 

0 

-1 

0 

-ei 

0 

e3 

ei3 

e23 

0 

1 

0 

0 

~e2 

-e3 

0 

e23 

ei23 

1 

ei 

e2 

e3 

ei2 

ei3 

e23 

ei23 
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Appendix  B  Listing  of  Clifford  algebra. m 

7.CA  =  CLIFFORD_ALGEBRA  (n) 

7. 

^Clifford  algebra  of  an  n-dimensional  vector  space.  Hultivectors  are 
70represented  by  2~n-column  matrices.  All  algebraic  operations  are 
70vectorized. 

7. 

7,  d  =  CA.  dimension  -  algebra  dimension  (d  =  2~n) 

7,  c  =  CA. set_scalar(s)  -  create  multivector  from  scalar 

7,  s  =  CA.get_scalar(c)  -  extract  scalar  from  multivector 

7,  c  =  CA. set_vector(v)  -  create  multivector  from  vector 

7,  v  =  CA.get_vector(c)  -  extract  vector  from  multivector 

7,  c  =  CA.even(c)  -  projection  to  the  even  subalgebra 

7,  c  =  CA. reverse (c)  -  reversion  operator  on  multivector 

7.  c  =  CA.dual(c)  -  dual  operator  on  multivector 

7,  c  =  CA.dual_inv(c)  -  inverse  dual  operator  on  multivector 

7.  c  =  CA. product  (a, b)  -  geometric  product  of  two  multivectors 

7,  c  =  CA.product_p(a,b)  -  progressive  outer  product  of  two  multivectors 

7,  c  =  CA.product_r (a,b)  -  regressive  outer  product  of  two  multivectors 

7.  c  =  CA.product_s(a,b)  -  inner  (scalar)  product  of  two  multivectors 

7,  CA.displayO  -  display  basic  information 

7. 

7#Copyright  (C)  2014  Defence  Science  and  Technology  Organisation 
7. 

7oCreated  by  Leonid  K.  Antanovskii 
function  CA  =  cliff ord_algebra(n) 
switch  n 

case  {1,2,3,4} 
otherwise 

error(’Space  dimension  7»d  not  supported.  5  ,n)  ; 

end 

d  =  2~n; 

function  check_dimension(c) 
if  size(c,2)  ~=  d 

error(’Wrong  multivector  dimension.’); 

end 

end 

function  c  =  set_scalar (s) 
if  size(s,2)  ~=  1 

error (’Wrong  scalar  dimension.’); 

end 

c  =  zeros (size (s , 1) ,d) ; 
c(:  ,1)  =  s; 

end 

function  c  =  set_vector (v) 
if  size(v,2)  ~=  n 

error (’Wrong  vector  dimension.’); 

end 

c  =  zeros (size (v, 1) ,d) ; 
c( : ,2 :n+l)  =  v; 

end 

function  s  =  get_scalar (c) 
check_dimension(c) ; 
s  =  c(: ,1) ; 

end 

function  v  =  get_vector (c) 
check_dimension(c) ; 
v  =  c( : ,2 :n+l) ; 

end 
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function  c  =  reverse (c) 
check_dimension(c) ; 
switch  n 
case  1 
case  2 

c  =  [c ( : ,1:3) ,-c(: ,4)] ; 
case  3 

c  =  [c(: ,1:4) ,-c(: ,5:8)] ; 
case  4 

c  =  [c(: ,1:5) ,-c(: ,6:15) ,c(: ,16)] ; 

end 

end 

function  c  =  even(c) 

check_dimension(c) ; 
switch  n 
case  1 

c( : ,2)  =  0.0; 
case  2 

c ( : ,  [2,3])  =  0.0; 
case  3 

c ( : ,  [2, 3, 4, 8])  =  0.0; 

case  4 

c(: , [2,3,4,5,12,13,14,15])  =  0.0; 

end 

end 

function  c  =  dual(c) 

check_dimension(c) ; 
switch  n 
case  1 

c  =  dualld(c)  ; 
case  2 

c  =  dual2d(c)  ; 
case  3 

c  =  dual3d(c)  ; 
case  4 

c  =  dual4d(c)  ; 

end 

end 

function  c  =  dual_inv(c) 
check_dimension(c) ; 
switch  n 
case  1 

c  =  dualld(c)  ; 
case  2 

c  =  dual2d_inv(c) ; 
case  3 

c  =  dual3d(c)  ; 
case  4 

c  =  dual4d_inv(c) ; 

end 

end 

function  c  =  product (a, b) 
check_dimension(a) ; 
check_dimension(b) ; 
switch  n 
case  1 

c  =  product Id (a, b) ; 
case  2 

c  =  product2d(a,b) ; 
case  3 

c  =  product3d(a,b) ; 
case  4 

c  =  product4d(a,b) ; 
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end 

end 

function  c  =  product_p(a,b) 
check_dimension(a) ; 
check_dimension(b) ; 
switch  n 
case  1 

c  =  product ld_p (a, b) ; 
case  2 

c  =  product2d_p(a,b) ; 
case  3 

c  =  product3d_p(a,b) ; 
case  4 

c  =  product4d_p(a,b) ; 

end 

end 

function  c  =  product_r (a,b) 
check_dimension(a) ; 
check_dimension(b) ; 
switch  n 
case  1 

c  =  productld_r (a,b) ; 
case  2 

c  =  product2d_r (a,b) ; 
case  3 

c  =  product3d_r (a,b) ; 
case  4 

c  =  product4d_r (a,b) ; 

end 

end 

function  c  =  product_s (a,b) 
check_dimension(a) ; 
check_dimension(b) ; 
switch  n 
case  1 

c  =  productld(a,b) ; 
case  2 

c  =  product2d_s (a,b) ; 
case  3 

c  =  product3d_s (a,b) ; 
case  4 

c  =  product4d_s (a,b) ; 

end 

end 


function  display 

fprintf  ( ’Function:  #/*s\n’  ,mf  ilename)  ; 
fprintf  (  ’Space  dimension:  °/0d\n’ ,n)  ; 
fprintf  (  ’Algebra  dimension:  0/0d\n’,d); 
fprintf ( ’Basis  multivectors : \n’ ) ; 
switch  n 
case  1 

fprintf ( ’\te(l)  =1  -  scalar\n’); 


fprintf ( ’\te(2) 

!  2 

=  el 

-  vector  (or  pseudoscalar) \n: 

fprintf ( ’\te(l) 

=  1 

-  scalar\n’); 

fprintf ( ’\te(2) 

=  el 

-  vector  l\n’); 

fprintf ( ’ \te (3) 

=  e2 

-  vector  2\n’); 

fprintf ( ’\te(4) 

!  3 

=  el2 

-  pseudoscalar \n’ ) ; 

fprintf ( ’\te(l) 

=  1 

-  scalar\n’); 

fprintf ( ’\te(2) 

-  el 

-  vector  l\n’); 

fprintf ( ’ \te (3) 

=  e2 

-  vector  2\n’); 

fprintf ( ’\te(4) 

=  e3 

-  vector  3\n’); 
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fprintf ( ’ \te(5) 
fprintf ( ’XteCe) 
fprintf ( ’XteC?) 
fprintf ( ’ \te(8) 
case  4 

fprintf ( ’YteCl) 
fprintf ( ,\te(2) 
fprintf ( ’ \te (3) 
fprintf ( ,\te(4) 
fprintf ( ’XteCS) 
fprintf ( ’ \te (6) 
fprintf ( ’ \te(7) 
fprintf ( ’ \te(8) 
fprintf ( *\te(9) 
fprintf ( ’XteClO) 
fprintf ( 5\te(ll) 
fprintf ( ,\te(12) 
fprintf ( ’ \te(13) 
fprintf ( J\te(14) 
fprintf ( ’ \te(15) 
fprintf ( ’ \te (16) 

end 

end 

CA. dimension  =  d; 

CA. set_scalar  =  ©set_scalar; 

CA.get_scalar  =  @get_scalar ; 

CA. set_vector  =  @set_vector; 

CA.get_vector  =  @get_vector; 

CA.even  =  ©even; 

CA. reverse  =  ©reverse; 

CA.dual  =  ©dual; 

CA.dual_inv  =  ©dual_inv; 

CA. product  =  ©product; 

CA.product_p  =  @product_p; 

CA.product_r  =  ©product_r; 

CA.product_s  =  ©product_s; 

CA. display  =  ©display; 

end 


=  el2 

bi-vector  l\nJ); 

=  el3 

bi-vector  2\n’); 

=  e23 

bi-vector  3\n’); 

=  el23  - 

pseudoscalar \n’ ) ; 

=  1 

-  scalar\nJ); 

=  el 

-  vector  l\n’ ) ; 

=  e2 

-  vector  2\n’ ) ; 

=  e3 

-  vector  3\n’ ) ; 

=  e4 

-  vector  4\n’ ) ; 

=  el2 

-  bi-vector  l\n’) 

=  el3 

-  bi-vector  2\n’) 

=  el4 

-  bi-vector  3\n’) 

=  e23 

-  bi-vector  4\n’) 

=  e24 

-  bi-vector  5\n’) 

=  e34 

-  bi-vector  6\n’) 

=  el23 

-  tri-vector  l\n’) 

=  el24 

-  tri-vector  2\n’) 

=  el34 

-  tri-vector  3\n’) 

=  e234 

-  tri-vector  4\n’) 

=  el234 

-  pseudoscalarXn’ ) 

%===================== 

function  c  =  dualld(c) 
cO  =  c( : , 1) ; 
cl  =  c( : ,2) ; 
c  =  [cl , cO] ; 
end 


7«========================== 

function  c  =  product Id (a ,b) 

aO  =  a( : , 1) ; 

al  =  a( : ,2) ; 

bO  =  b(: ,1) ; 

bl  =  b(: ,2) ; 

cO  =  aO.*bO  +  al.*bl; 

cl  =  aO.*bl  +  al.*bO; 

c  =  [cO,cl] ; 

end 


7.============================ 

function  c  =  product ld_p (a, b) 
aO  =  a( : , 1) ; 
al  =  a(: ,2) ; 

bO  =  b(: ,1) ; 
bl  =  b(: ,2) ; 

cO  =  aO.*bO; 
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cl  =  aO.*bl  +  al.*bO; 

c  =  [cO,cl] ; 

end 

%============================ 

function  c  =  product ld_r (a, b) 

a  =  dualld(a) ; 

b  =  dual Id (b) ; 

c  =  productld_p(a,b) ; 

c  =  dual Id (c) ; 

end 

%============================ 


function 

c  - 

dual 2d (c) 

cO  = 

c( : 

i); 

cl  = 

c( : 

2); 

c2  = 

c( : 

3); 

cl2 

=  c( 

,4) 

c  = 

[cl2 

c2. 

-cl , cO] ; 

end 

l- 


function 

c  =  dual2d_inv(c) 

cO  = 

c( : 

i); 

cl  = 

c( : 

2); 

c2  = 

c( : 

3); 

cl2 

=  c( 

,4); 

c  = 

[cl2 

-c2, cl , cO] ; 

end 

%===== - ================== - ====== 

function  c  =  product 2d (a, b) 
aO  =  a( : , 1) ; 
al  =  a( : ,2) ; 
a2  =  a( : ,3) ; 
al2  =  a(: ,4); 

bO  =  b(: ,1) ; 
bl  =  b(: ,2) ; 
b2  =  b(: ,3) ; 
bl2  =  b(: ,4); 

cO  =  aO . *b0  +  al . *bl  +  a2.*b2  -  al2.*bl2; 
cl  =  aO.*bl  +  al.*b0  -  a2.*bl2  +  al2.*b2; 
c2  =  a0.*b2  +  al.*bl2  +  a2.*b0  -  al2.*bl; 
cl2  =  a0.*bl2  +  al.*b2  -  a2.*bl  +  al2.*b0; 

c  =  [cO,cl ,c2,cl2] ; 
end 

%========================================= 

function  c  =  product2d_p(a,b) 
aO  =  a( : , 1) ; 
al  =  a( : ,2) ; 
a2  =  a( : ,3) ; 
al2  =  a(: ,4); 

bO  =  b(: ,1) ; 
bl  =  b(: ,2) ; 
b2  =  b(: ,3) ; 
bl2  =  b(: ,4); 

cO  =  aO.*bO; 

cl  =  aO.*bl  +  al.*bO; 

c2  =  a0.*b2  +  a2.*b0; 

cl2  =  a0.*bl2  +  al.*b2  -  a2.*bl  +  al2.*b0; 
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c  =  [cO,cl ,c2,cl2] ; 
end 


%============================ 

function  c  =  product2d_r (a,b) 

a  =  dual2d(a) ; 

b  =  dual2d(b) ; 

c  =  product2d_p(a,b) ; 

c  =  dual2d_inv(c) ; 

end 


%=======================================: 

function  c  =  product2d_s (a,b) 
aO  =  a( : , 1) ; 
al  =  a( : ,2) ; 
a2  =  a( : ,3) ; 
al2  =  a(: ,4); 

bO  =  b(: ,1) ; 
bl  =  b(: ,2) ; 
b2  =  b(: ,3) ; 
bl2  =  b(: ,4) ; 

cO  =  aO.*bO  +  al.*bl  +  a2.*b2  -  al2.*bl2 
cl  =  aO.*bl  +  al.*bO  -  a2.*bl2  +  al2.*b2 
c2  =  a0.*b2  +  al.*bl2  +  a2.*b0  -  al2.*bl 
cl2  =  aO . *bl2  +  al2 . *b0; 

c  =  [cO,cl ,c2,cl2] ; 
end 


l- 


function 

c  =  dual3d(c) 

cO  = 

c( : 

1); 

cl  = 

c( : 

2); 

c2  = 

c( : 

3); 

c3  = 

c  ( 

4); 

cl2 

=  c( 

,5); 

cl3 

=  c( 

,6); 

c23 

=  c( 

,7); 

cl23 

=  c 

:  ,8)  ; 

c  =  [cl23 , c23 , -cl3 , cl2 , c3 , -c2 , cl , cO] ; 
end 


7o========================== 

function  c  =  product3d(a,b) 


aO  = 

a( : 

i); 

al  = 

a( : 

2); 

a2  = 

a( : 

3); 

a3  = 

a( : 

4); 

al2 

=  a( 

,5); 

al3 

=  a( 

,6); 

a23 

=  a( 

,7); 

al23 

=  a( 

:  ,8) 

bO  =  b(: ,1) ; 
bl  =  b(: ,2) ; 
b2  =  b(: ,3) ; 
b3  =  b(: ,4) ; 
bl2  =  b( : ,5) ; 
bl3  =  b( : ,6) ; 
b23  =  b(: ,7) ; 
bl23  =  b( : ,8) ; 

cO  =  aO . *b0  . . . 

+  al.*bl  +  a2.*b2  +  a3.*b3  ... 

-  al2 . *bl2  -  al3 . *bl3  -  a23.*b23  ... 
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-  al23.*bl23; 
cl  =  aO . *bl  . . . 

+  al.*b0  -  a2.*bl2  -  a3.*bl3  ... 

+  al2 .  *b2  +  al3.*b3  -  a23.*bl23  ... 

-  al23.*b23; 
c2  =  aO . *b2  .  .  . 

+  al.*bl2  +  a2.*b0  -  a3.*b23  ... 

-  al2.*bl  +  al3.*bl23  +  a23.*b3  ... 
+  al23 . *bl3; 

c3  =  aO . *b3  .  .  . 

+  al.*bl3  +  a2.*b23  +  a3.*b0  ... 

-  al2 . *bl23  -  al3.*bl  -  a23.*b2  ... 

-  al23.*bl2; 
cl2  =  aO . *bl2  . . . 

+  al.*b2  -  a2.*bl  +  a3.*bl23  ... 

+  al2 .  *b0  -  al3.*b23  +  a23.*bl3  ... 
+  al23.*b3; 
cl3  =  aO . *bl3  . .  . 

+  al . *b3  -  a2.*bl23  -  a3.*bl  ... 

+  al2 . *b23  +  al3 . *bO  -  a23.*bl2  ... 

-  al23 . *b2 ; 
c23  =  aO . *b23  . . . 

+  al.*bl23  +  a2.*b3  -  a3.*b2  ... 

-  al2 . *bl3  +  al3 . *bl2  +  a23.*b0  ... 
+  al23.*bl; 

cl23  =  aO . *bl23  . .  . 

+  al . *b23  -  a2.*bl3  +  a3.*bl2  ... 

+  al2 . *b3  -  al3 . *b2  +  a23.*bl  ... 

+  al23.*b0; 

c  =  [00,01,02,03,012,013,023,0123]; 
end 


/o 

function 

c  =  product3d_p(a,b) 

aO  = 

a( : 

l); 

al  = 

a( : 

2); 

a2  = 

a( : 

3); 

a3  = 

a( : 

4); 

al2 

=  a( 

,5); 

al3 

=  a( 

,6); 

a23 

=  a( 

,7); 

al23 

=  a( 

:  ,8)  ; 

bO  =  b(: ,1) ; 
bl  =  b(: ,2) ; 
b2  =  b(: ,3) ; 
b3  =  b(: ,4) ; 
bl2  =  b(: ,5) ; 
bl3  =  b( : ,6) ; 
b23  =  b(: ,7) ; 
bl23  =  b(: ,8) ; 

cO  =  aO.*bO; 
cl  =  aO.*bl  +  al.*bO; 
c2  =  a0.*b2  +  a2.*b0; 
c3  =  a0.*b3  +  a3.*b0; 

cl2  =  a0.*bl2  +  al.*b2  -  a2.*bl  +  al2.*b0; 

cl3  =  a0.*bl3  +  al.*b3  -  a3.*bl  +  al3.*b0; 

c23  =  aO . *b23  +  a2.*b3  -  a3.*b2  +  a23.*b0; 

cl23  =  aO . *bl23  . . . 

+  al . *b23  -  a2.*bl3  +  a3.*bl2  ... 

+  al2 . *b3  -  al3.*b2  +  a23.*bl  ... 

+  al23.*b0; 

c  =  [cO, cl , c2 , c3, cl2 , cl3, c23 , cl23] ; 
end 
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•/,============================ 

function  c  =  product3d_r (a,b) 

a  =  dual3d(a); 

b  =  dual3d(b) ; 

c  =  product3d_p(a,b) ; 

c  =  dual3d(c); 

end 


•/,============================ 

function  c  =  product3d_s(a,b) 


aO  = 

a( : 

i); 

al  = 

a( : 

2); 

a2  = 

a( : 

3); 

a3  = 

a( : 

4); 

al2 

=  a( 

,5); 

al3 

=  a( 

,6); 

a23 

=  a( 

,7); 

al23 

=  a( 

:  ,8) 

bO  =  b(: ,1) ; 
bl  =  b( : ,2) ; 
b2  =  b(: ,3) ; 
b3  =  b(: ,4) ; 
bl2  =  b( : ,5) ; 
bl3  =  b( : ,6) ; 
b23  =  b(: ,7) ; 
bl23  =  b(: ,8) ; 

cO  =  aO.*bO  +  al.*bl  +  a2.*b2  +  a3.*b3  ... 

-  al2 . *bl2  -  al3 . *bl3  -  a23.*b23  -  al23.*bl23; 
cl  =  aO.*bl  +  al.*bO  -  a2.*bl2  -  a3.*bl3  ... 

+  al2 . *b2  +  al3.*b3  -  a23.*bl23  -  al23.*b23; 
c2  =  a0.*b2  +  al.*bl2  +  a2.*b0  -  a3.*b23  ... 

-  al2 . *bl  +  al3 . *bl23  +  a23.*b3  +  al23.*bl3; 
c3  =  a0.*b3  +  al.*bl3  +  a2.*b23  +  a3.*b0  ... 

-  al2 . *bl23  -  al3.*bl  -  a23.*b2  -  al23.*bl2; 

cl2  =  a0.*bl2  +  a3 . *bl23  +  al2 . *b0  +  al23.*b3; 

cl3  =  aO .  *bl3  -  a2 .  *bl23  +  al3.*b0  -  al23.*b2; 

c23  =  aO .  *b23  +  al.*bl23  +  a23.*b0  +  al23.*bl; 

cl23  =  aO . *bl23  +  al23.*b0; 

c  =  [cO, cl , c2 , c3, cl2 , cl3, c23 , cl23] ; 
end 


•/,===================== 

function  c  =  dua!4d(c) 


cO 

= 

c( 

,1) 

cl 

= 

c( 

,2) 

c2 

= 

c( 

,3) 

c3 

= 

c( 

,4) 

c4 

= 

c( 

,5) 

cl2  =  c ( : ,6) ; 
cl3  =  c(:,7); 
cl4  =  c(: ,8) ; 
c23  =  c ( : ,9) ; 
c24  =  c(: ,10); 
c34  =  c(: ,11) ; 
cl23  =  c(: ,12) ; 
cl24  =  c(: ,13) ; 
cl34  =  c(: ,14) ; 
c234  =  c(: ,15) ; 
cl 234  =  c( : ,16); 
c  =  [cl234, . . . 

c234,-cl34,cl24,-cl23, . . . 
c34 , -c24 , c23 , cl4 , -cl3 , cl2 , . . . 
c4,-c3, c2 , -cl , . . . 
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c0]  ; 

end 


function 

c  =  dual4d_inv(c) 

cO  =  c( : 

i); 

cl  =  c( : 

2); 

c2  =  c( : 

3); 

c3  =  c( : 

4); 

c4  =  c( : 

5); 

cl2  =  c ( 

,6); 

cl3  =  c ( 

,7  ); 

cl4  =  c( 

,8); 

c23  =  c ( 

,9); 

c24  =  c ( 

,10); 

c34  =  c ( 

,n); 

cl23  =  c(: ,12) ; 

cl24  =  c(: ,13) ; 

cl34  =  c(: ,14) ; 

c234  =  c(: ,15) ; 

cl 234  =  c(: ,16) ; 

c  =  [cl234, . . . 

-c234,cl34,-cl24,cl23, . . . 
034,-024,023,014,-013,012, . . . 
-c4, c3, -c2 , cl , .  .  . 
cO]  ; 

end 


•/,======================================================================== 

function  c  =  product4d(a,b) 

aO  =  a(  :  ,  1)  ; 

al  =  a(  :  ,2)  ; 

a2  =  a(: ,3) ; 

a3  =  a(  :  ,4)  ; 

a4  =  a(  :  ,5)  ; 

al2  =  a( : ,6) ; 

al3  =  a( : ,7) ; 

al4  =  a( : ,8) ; 

a23  =  a( : , 9) ; 

a24  =  a(: ,10); 

a34  =  a(: ,11); 

al23  =  a(: ,12) ; 

al24  =  a(: ,13) ; 

al34  =  a(: ,14) ; 

a234  =  a(: ,15) ; 

al234  =  a(: ,16); 

bO  =  b( : ,1) ; 
bl  =  b( : ,2) ; 
b2  =  b(: ,3) ; 
b3  =  b(: ,4) ; 
b4  =  b(: ,5) ; 
bl2  =  b( : ,6) ; 
bl3  =  b(: ,7) ; 
bl4  =  b( : ,8) ; 
b23  =  b( : ,9) ; 
b24  =  b(: ,10); 
b34  =  b(: ,11) ; 
bl23  =  b(: ,12) ; 
bl24  =  b(: ,13) ; 
bl34  =  b(: ,14) ; 
b234  =  b(: ,15) ; 
bl234  =  b( : ,16); 

cO  =  aO . *bO  . . . 

+  al.*bl  +  a2.*b2  +  a3.*b3  +  a4.*b4  ... 

-  al2 . *bl2  -  al3 . *bl3  -  al4.*bl4  -  a23.*b23  -  a24.*b24  -  a34.*b34  ... 
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-  al23 .  *bl23  -  al24.*bl24  -  al34.*bl34  -  a234.*b234  ... 

+  al234 . *bl234 ; 

cl  =  aO.*bl  +  al.*bO  -  a2.*bl2  -  a3.*bl3  -  a4.*bl4  ... 

+  al2 .  *b2  +  al3.*b3  +  al4.*b4  -  a23 . *bl23  -  a24.*bl24  -  a34.*bl34 

-  al23.*b23  -  al24.*b24  -  al34.*b34  +  a234.*bl234  ... 

-  al234.*b234; 
c2  =  aO . *b2  . . . 

+  al.*bl2  +  a2.*b0  -  a3.*b23  -  a4.*b24  ... 

-  al2 . *bl  +  al3.*bl23  +  al4.*bl24  +  a23.*b3  +  a24.*b4  -  a34.*b234 
+  al23.*bl3  +  al24.*bl4  -  al34.*bl234  -  a234.*b34  ... 

+  al234.*bl34; 
c3  =  aO . *b3  .  .  . 

+  al.*bl3  +  a2.*b23  +  a3.*b0  -  a4.*b34  ... 

-  al2 . *bl23  -  al3 . *bl  +  al4.*bl34  -  a23.*b2  +  a24.*b234  +  a34.*b4 

-  al23.*bl2  +  al24.*bl234  +  al34.*bl4  +  a234.*b24  ... 

-  al234 . *bl24; 
c4  =  aO . *b4  . . . 

+  al.*bl4  +  a2.*b24  +  a3.*b34  +  a4.*b0  ... 

-  al2 . *bl24  -  al3 . *bl34  -  al4.*bl  -  a23.*b234  -  a24.*b2  -  a34.*b3 

-  al23.*bl234  -  al24.*bl2  -  al34.*bl3  -  a234.*b23  ... 

+  al234 .  *bl23 ; 

cl2  =  aO . *bl2  . . . 

+  al.*b2  -  a2.*bl  +  a3.*bl23  +  a4.*bl24  ... 

+  al2 .  *bO  -  al3.*b23  -  al4.*b24  +  a23.*bl3  +  a24.*bl4  -  a34.*bl234 
+  al23 .  *b3  +  al24 .  *b4  -  al34.*b234  +  a234.*bl34  ... 

-  al234 .  *b34; 
cl3  =  aO . *bl3  . . . 

+  al.*b3  -  a2.*bl23  -  a3.*bl  +  a4.*bl34  ... 

+  al2 . *b23  +  al3 . *bO  -  al4.*b34  -  a23.*bl2  +  a24.*bl234  +  a34.*bl4 

-  al23 . *b2  +  al24 . *b234  +  al34.*b4  -  a234.*bl24  ... 

+  al234.*b24; 

cl4  =  aO . *bl4  . . . 

+  al . *b4  -  a2.*bl24  -  a3.*bl34  -  a4.*bl  ... 

+  al2 . *b24  +  al3 . *b34  +  al4.*b0  -  a23.*bl234  -  a24.*bl2  -  a34.*bl3 

-  al23.*b234  -  al24.*b2  -  al34.*b3  +  a234.*bl23  ... 

-  al234.*b23; 
c23  =  aO . *b23  . . . 

+  al.*bl23  +  a2.*b3  -  a3.*b2  +  a4.*b234  ... 

-  al2 . *bl3  +  al3 . *bl2  -  al4.*bl234  +  a23.*b0  -  a24.*b34  +  a34.*b24 

+  al23.*bl  -  al24 . *bl34  +  al34.*bl24  +  a234.*b4  ... 

-  al234.*bl4; 
c24  =  aO . *b24  . . . 

+  al . *bl24  +  a2 . *b4  -  a3.*b234  -  a4.*b2  ... 

-  al2 . *bl4  +  al3 . *bl234  +  al4.*bl2  +  a23.*b34  +  a24.*b0  -  a34.*b23 

+  al23.*bl34  +  al24.*bl  -  al34.*bl23  -  a234.*b3  ... 

+  al234 .  *bl3; 
c34  =  aO . *b34  . . . 

+  al .  *bl34  +  a2.*b234  +  a3.*b4  -  a4.*b3  ... 

-  al2 .  *bl234  -  al3.*bl4  +  al4.*bl3  -  a23.*b24  +  a24.*b23  +  a34.*b0 

-  al23.*bl24  +  al24.*bl23  +  al34.*bl  +  a234.*b2  ... 

-  al234.*bl2; 
cl23  =  aO . *bl23  . . . 

+  al .  *b23  -  a2.*bl3  +  a3.*bl2  -  a4.*bl234  ... 

+  al2 . *b3  -  al3.*b2  +  al4.*b234  +  a23.*bl  -  a24.*bl34  +  a34.*bl24 

+  al23.*b0  -  al24.*b34  +  al34.*b24  -  a234.*bl4  ... 

+  al234 .  *b4 ; 
cl24  =  aO . *bl24  . . . 

+  al.*b24  -  a2 . *bl4  +  a3.*bl234  +  a4.*bl2  ... 

+  al2 . *b4  -  al3.*b234  -  al4.*b2  +  a23.*bl34  +  a24.*bl  -  a34.*bl23 

+  al23.*b34  +  al24.*b0  -  al34.*b23  +  a234.*bl3  ... 

-  al234.*b3; 
cl34  =  aO . *bl34  . . . 

+  al . *b34  -  a2.*bl234  -  a3.*bl4  +  a4.*bl3  ... 

+  al2 . *b234  +  al3.*b4  -  al4.*b3  -  a23.*bl24  +  a24.*bl23  +  a34.*bl 

-  al23.*b24  +  al24.*b23  +  al34.*b0  -  a234.*bl2  ... 

+  al234 . *b2 ; 

c234  =  aO . *b234  . . . 
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+  al.*bl234  +  a2.*b34  -  a3.*b24  +  a4.*b23  ... 

-  al2 . *bl34  +  al3.*bl24  -  al4.*bl23  +  a23.*b4  -  a24.*b3  +  a34.*b2  ... 

+  al23 . *bl4  -  al24.*bl3  +  al34.*bl2  +  a234.*b0  ... 

-  al234.*bl; 

cl 234  =  aO . *bl234  . . . 

+  al .  *b234  -  a2.*bl34  +  a3.*bl24  -  a4.*bl23  ... 

+  al2 . *b34  -  al3 . *b24  +  al4.*b23  +  a23.*bl4  -  a24.*bl3  +  a34.*bl2  ... 

+  al23.*b4  -  al24 . *b3  +  al34.*b2  -  a234.*bl  +  al234.*b0; 

c  =  [c0,cl,c2,c3,c4,cl2,cl3,cl4,c23,c24,c34,cl23,cl24,cl34,c234,cl234] ; 
end 

•/,======================================================================== 

function  c  =  product4d_p(a,b) 

aO  =  a(  :  ,  1)  ; 

al  =  a(  :  ,2)  ; 

a2  =  a( : ,3) ; 

a3  =  a(  :  ,4)  ; 

a4  =  a(: ,5) ; 

al2  =  a( : ,6) ; 

al3  =  a( : ,7) ; 

al4  =  a( : ,8) ; 

a23  =  a( : ,9) ; 

a24  =  a(: ,10); 

a34  =  a(: ,11); 

al23  =  a(: ,12) ; 

al24  =  a(: ,13) ; 

al34  =  a(: ,14) ; 

a234  =  a(: ,15) ; 

al234  =  a(: ,16); 

bO  =  b( : ,1) ; 
bl  =  b(: ,2) ; 
b2  =  b(: ,3) ; 
b3  =  b(: ,4) ; 
b4  =  b(: ,5) ; 
bl2  =  b( : ,6) ; 
bl3  =  b(: ,7) ; 
bl4  =  b( : ,8) ; 
b23  =  b( : ,9) ; 
b24  =  b(: ,10); 
b34  =  b(: ,11); 
bl23  =  b(: ,12) ; 
bl24  =  b(: ,13) ; 
bl34  =  b(: ,14) ; 
b234  =  b(: ,15) ; 
bl234  =  b(: ,16) ; 

cO  =  aO.*bO; 
cl  =  aO.*bl  +  al.*bO; 
c2  =  a0.*b2  +  a2.*b0; 
c3  =  a0.*b3  +  a3.*b0; 
c4  =  aO .  *b4  +  a4 .  *b0 ; 

cl2  =  a0.*bl2  +  al.*b2  -  a2.*bl  +  al2.*b0; 

cl3  =  a0.*bl3  +  al.*b3  -  a3.*bl  +  al3.*b0; 

cl4  =  a0.*bl4  +  al.*b4  -  a4.*bl  +  al4.*b0; 

c23  =  a0.*b23  +  a2.*b3  -  a3.*b2  +  a23.*b0; 

c24  =  a0.*b24  +  a2.*b4  -  a4.*b2  +  a24.*b0; 

c34  =  aO . *b34  +  a3.*b4  -  a4.*b3  +  a34.*b0; 

cl23  =  aO . *bl23  +  al . *b23  -  a2.*bl3  +  a3.*bl2  ... 

+  al2 . *b3  -  al3.*b2  +  a23.*bl  +  al23.*b0; 
cl24  =  aO . *bl24  +  al . *b24  -  a2.*bl4  +  a4.*bl2  ... 

+  al2 . *b4  -  al4 . *b2  +  a24.*bl  +  al24.*b0; 
cl34  =  aO .  *bl34  +  al .  *b34  -  a3.*bl4  +  a4.*bl3  ... 

+  al3.*b4  -  al4.*b3  +  a34.*bl  +  al34.*b0; 
c234  =  aO . *b234  +  a2.*b34  -  a3.*b24  +  a4.*b23  ... 

+  a23.*b4  -  a24.*b3  +  a34.*b2  +  a234.*b0; 
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cl 234  =  a0.*bl234  .  .  . 

+  al.*b234  -  a2.*bl34  +  a3.*bl24  -  a4.*bl23  ... 

+  al2 . *b34  -  al3 . *b24  +  al4.*b23  +  a23.*bl4  -  a24.*bl3  +  a34.*bl2  ... 

+  al23.*b4  -  al24 . *b3  +  al34.*b2  -  a234.*bl  +  al234.*b0; 

c  =  [c0,cl,c2,c3,c4,cl2,cl3,cl4,c23,c24,c34,cl23,cl24,cl34,c234,cl234] ; 
end 

•/,======================================================================== 

function  c  =  product4d_r (a,b) 

a  =  dual4d(a); 

b  =  dual4d(b); 

c  =  product4d_p(a,b) ; 

c  =  dual4d_inv(c) ; 

end 

•/,======================================================================== 

function  c  =  product4d_s(a,b) 

aO  =  a(  :  ,  1)  ; 

al  =  a(  :  ,2)  ; 

a2  =  a( : ,3) ; 

a3  =  a(  :  ,4)  ; 

a4  =  a(  :  ,5)  ; 

al2  =  a( : , 6) ; 

al3  =  a(:,7); 

al4  =  a( : ,8) ; 

a23  =  a( : ,9) ; 

a24  =  a(: ,10); 

a34  =  a(: ,11) ; 

al23  =  a(: ,12) ; 

al24  =  a(: ,13) ; 

al34  =  a(: ,14) ; 

a234  =  a(: ,15) ; 

al234  =  a(: ,16); 

bO  =  b( : ,1) ; 
bl  =  b( : ,2) ; 
b2  =  b(: ,3) ; 
b3  =  b(: ,4) ; 
b4  =  b(: ,5) ; 
bl2  =  b( : ,6) ; 
bl3  =  b(: ,7) ; 
bl4  =  b( : ,8) ; 
b23  =  b( : ,9) ; 
b24  =  b(: ,10); 
b34  =  b(: ,11); 
bl23  =  b(: ,12) ; 
bl24  =  b(: ,13) ; 
bl34  =  b(: ,14) ; 
b234  =  b(: ,15) ; 
bl234  =  b( : ,16); 

cO  =  aO . *b0  . . . 

+  al.*bl  +  a2.*b2  +  a3.*b3  +  a4.*b4  ... 

-  al2 . *bl2  -  al3 . *bl3  -  al4.*bl4  -  a23.*b23  -  a24.*b24  -  a34.*b34  ... 

-  al23.*bl23  -  al24.*bl24  -  al34.*bl34  -  a234.*b234  ... 

+  al234.*bl234; 

cl  =  aO.*bl  +  al.*b0  -  a2.*bl2  -  a3.*bl3  -  a4.*bl4  ... 

+  al2 . *b2  +  al3.*b3  +  al4.*b4  -  a23 . *bl23  -  a24.*bl24  -  a34.*bl34  ... 

-  al23.*b23  -  al24.*b24  -  al34.*b34  +  a234.*bl234  ... 

-  al234.*b234; 
c2  =  aO . *b2  .  .  . 

+  al.*bl2  +  a2.*b0  -  a3.*b23  -  a4.*b24  ... 

-  al2 . *bl  +  al3.*bl23  +  al4.*bl24  +  a23.*b3  +  a24.*b4  -  a34.*b234  ... 
+  al23 .  *bl3  +  al24.*bl4  -  al34.*bl234  -  a234.*b34  ... 

+  al234.*bl34; 
c3  =  aO . *b3  . . . 
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+  al.*bl3  +  a2.*b23  +  a3.*b0  -  a4.*b34  ... 

-  al2 . *bl23  -  al3.*bl  +  al4.*bl34  -  a23.*b2  +  a24.*b234  +  a34.*b4  ... 

-  al23 . *bl2  +  al24.*bl234  +  al34.*bl4  +  a234.*b24  ... 

-  al234 .  *bl24; 
c4  =  aO . *b4  . . . 

+  al.*bl4  +  a2.*b24  +  a3.*b34  +  a4.*b0  ... 

-  al2 . *bl24  -  al3 . *bl34  -  al4.*bl  -  a23.*b234  -  a24.*b2  -  a34.*b3  ... 

-  al23.*bl234  -  al24.*bl2  -  al34.*bl3  -  a234.*b23  ... 

+  al234 .  *bl23 ; 

cl2  =  a0.*bl2  +  a3 . *bl23  +  a4.*bl24  +  al2 . *b0  -  a34.*bl234  ... 

+  al23 . *b3  +  al24 . *b4  -  al234.*b34; 
cl3  =  a0.*bl3  -  a2 . *bl23  +  a4.*bl34  +  al3.*b0  +  a24.*bl234  ... 

-  al23.*b2  +  al34.*b4  +  al234.*b24; 

cl4  =  a0.*bl4  -  a2 . *bl24  -  a3.*bl34  +  al4.*b0  -  a23.*bl234  ... 

-  al24.*b2  -  al34 . *b3  -  al234.*b23; 

c23  =  a0.*b23  +  al.*bl23  +  a4.*b234  -  al4.*bl234  +  a23.*b0  ... 

+  al23.*bl  +  a234.*b4  -  al234.*bl4; 
c24  =  a0.*b24  +  al.*bl24  -  a3.*b234  +  al3.*bl234  +  a24.*b0  ... 

+  al24 . *bl  -  a234.*b3  +  al234.*bl3; 
c34  =  a0.*b34  +  al.*bl34  +  a2.*b234  -  al2.*bl234  +  a34.*b0  ... 

+  al34 . *bl  +  a234.*b2  -  al234.*bl2; 
cl23  =  aO . *bl23  -  a4.*bl234  +  al23.*b0  +  al234.*b4; 

cl24  =  aO . *bl24  +  a3.*bl234  +  al24.*b0  -  al234.*b3; 

cl34  =  aO . *bl34  -  a2.*bl234  +  al34.*b0  +  al234.*b2; 

c234  =  aO .  *b234  +  al.*bl234  +  a234.*b0  -  al234.*bl; 

cl 234  =  aO .  *bl234  +  al234.*b0; 

c  =  [c0,cl,c2,c3,c4,cl2,cl3,cl4,c23,c24,c34,cl23,cl24,cl34,c234,cl234] ; 
end 
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Appendix  C  Listing  of  Clifford  algebra  test  .m 

‘/.CLIFFORD.  ALGEBRA_TEST  (N ,  tol) 

7. 

y.Unit  tests  of  basic  operations  of  Clifford  algebra. 

7. 

y.Copyright  (C)  2014  Defence  Science  and  Technology  Organisation 

7. 

y.Created  by  Leonid  K.  Antanovskii 
function  cliff ord_algebra_test (N,tol) 
fprintf ( ’Testing  Clifford  algebra\n’); 
if  nargin  <  2 

rng( ’default ’ ) ; 
tol  =  1.0e-14; 
if  nargin  <  1 
N  =  100000; 

end 

end 

for  n  =  1:4 

fprintf  ( ’Space  dimension:  7.d\n’ ,n)  ; 
test_geometric_product (n, tol ,N) ; 
test_progressive_product (n, tol ,N) ; 
test_regressive_product (n , tol , N) ; 
test_scalar_product (n,tol ,N) ; 
test_inner_product (n,tol ,N) ; 
test_even_subalgebra(n,N) ; 
test_reverse_operator (n,tol ,N) ; 
test_dual_operator (n,tol ,N) ; 
test_blade_f eature (n,tol ,N) ; 
test_identity_l (n,tol ,N) ; 
test_identity_2(n,tol ,N) ; 
test_identity_3(n,tol ,N) ; 
test_identity_4(n,tol ,N) ; 
test_identity_5(n,tol) ; 

end 

test_cross_product (tol ,N) ; 
test_complex_number (tol ,N) ; 
test_central_subalgebra(tol ,N) ; 
test_bicomplex_number (tol ,N) ; 
test_quaternion(tol,N) ; 
test.perf ormance(N) ; 

end 


•/,================================================================= 

function  test_geometric_product (n, tol ,N) 
fprintf ( ’Geometric  product  test\n’); 

CA  =  cliff ord.algebra(n) ; 
d  =  CA. dimension; 

7.  random  multivectors 
a  =  2.0*rand(N,d)  -  1.0; 
b  =  2.0*rand(N,d)  -  1.0; 
c  =  2.0*rand(N,d)  -  1.0; 

7.  check  associativity 

q  =  CA .product (CA .product (a, b) , c)  -  CA. product (a, CA. product (b,c)) ; 

assert (all (q( : )  <  tol)); 

end 


•/,=========================================== 

function  test_progressive_product (n,tol ,N) 
fprintf ( ’Progressive  outer  product  test\n’); 
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CA  =  cliff ord_algebra(n) ; 
d  =  CA. dimension; 

7.  random  multivectors 
a  =  2.0*rand(N,d)  -  1.0; 
b  =  2.0*rand(N,d)  -  1.0; 
c  =  2.0*rand(N,d)  -  1.0; 

7,  check  associativity 

q  =  CA.product_p(CA.product_p(a,b) ,c)  -  CA.product_p(a,CA.product_p(b,c)) ; 

assert (all (q( : )  <  tol)); 

end 


•/,========================================================================= 

function  test_regressive_product (n, tol ,N) 
fprintf ( ’Regressive  outer  product  test\n’); 

CA  =  cliff ord_algebra(n) ; 
d  =  CA. dimension; 

7.  random  multivectors 
a  =  2.0*rand(N,d)  -  1.0; 
b  =  2.0*rand(N,d)  -  1.0; 
c  =  2.0*rand(N,d)  -  1.0; 

7,  check  associativity 

q  =  CA.product_r(CA.product_r(a,b) ,c)  -  CA.product_r(a,CA.product_r(b,c)) ; 

assert (all (q( : )  <  tol)); 

end 


•/,==================================== 

function  test_scalar_product (n, tol ,N) 
fprintf ( ’Scalar  product  test\n’); 

CA  =  cliff ord_algebra(n) ; 

a  =  2 . 0*rand(N,n)  -  1.0; 
b  =  2 . 0*rand(N,n)  -  1.0; 
abl  =  CA.set_scalar(dot(a,b,2)) ; 

a  =  CA. set_vector (a) ; 
b  =  CA. set_vector (b) ; 
ab2  =  CA .product_s (a,b) ; 

assert (norm (abl  -  ab2,inf)  <  tol); 
end 


%================================================= 

function  test_inner_product (n , tol , N) 
fprintf (’ Inner  product  test\n’); 

CA  =  cliff ord_algebra(n) ; 

for  k  =  l:n 

fprintf  (’ \twedge  product  of  7»d  vectors\n’ ,k)  ; 

a  =  CA.set_vector(2.0*rand(N,n)  -  1.0); 
b  =  cell(k, 1) ; 
for  i  =  1 : k 

b{i}  =  CA . set_vector (2 . 0*rand(N,n)  -  1.0); 

end 

p  =  b{l>; 
for  i  =  2:k 

p  =  CA.product_p(p,b{i}) ; 

end 
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p  =  CA.product_s(a,p) ; 

q  =  CA. set_vector(zeros(N,n) ) ; 
for  i  =  1 : k 

s  =  CA.product_s(a,b{i}) ; 
for  j  =  l:k 
if  j  "=  i 

s  =  CA .product_p(s ,b{j}) ; 

end 

end 

q  =  q  +  (-1) ~ (i-1) *s ; 

end 

assert (norm (p  -  q,inf)  <  tol) ; 

end 

end 


•/======================================= 

function  test_even_subalgebra (n ,  N) 
fprintf ( ’Even  subalgebra  test\n’); 

CA  =  cliff ord_algebra(n) ; 
d  =  CA. dimension; 

7,  random  even  multivectors 
a  =  CA.even(2.0*rand(N,d)  -  1.0); 
b  =  CA.even(2.0*rand(N,d)  -  1.0); 

7*  check  subalgebra  property 
c  =  CA. product (a,b) ; 

assert (norm(CA. even(c)  -  c,inf)  ==  0.0); 
end 


%====================================== 

function  test_reverse_operator (n,tol ,N) 
fprintf ( ’Reversion  operator  test\n’); 

CA  =  cliff ord_algebra(n) ; 
d  =  CA. dimension; 

7.  random  multivectors 
a  =  2.0*rand(N,d)  -  1.0; 

7.  check  involution 

b  =  CA. reverse (a) ; 

c  =  CA. reverse (b) ; 

assert (norm(a  -  c,inf)  <  tol); 

end 


•/,=================================== 

function  test_dual_operator (n , tol , N) 
fprintf ( ’Dual  operator  test\n’); 

CA  =  cliff ord_algebra(n) ; 
d  =  CA. dimension; 

7.  random  multivectors 
a  =  2.0*rand(N,d)  -  1.0; 

7,  check  inversion 

b  =  CA.dual(a) ; 

c  =  CA.dual_inv(b) ; 

assert (norm(a  -  c,inf)  <  tol); 

end 


7o=================================== 

function  test_blade_f eature (n , tol , N) 
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fprintf ( ’Blade  feature  test\n’); 

CA  =  cliff ord_algebra(n) ; 
e  =  CA .  set_scalar  (1 . 0) ;  %  unit 

v  =  orth(2 . 0*rand(n,n)  -  1.0);  °/«  random  orthonormal  basis 
for  k  =  l:n 

f  printf  ( ’  \t°/0d-blade\n  ’ ,  k) ; 

°/0  construct  a  blade 

A  =  e; 

for  i  =  1 :k 

a  =  CA. set_vector (v(i , : ) ) ; 

A  =  CA. product (A, a) ; 

end 

7.  random  vectors 

a  =  CA. set_vector(2.0*rand(N,n)  -  1.0); 

%  check  blade  properties 

p  =  CA. product (a, A)  -  CA.product_s(a, A)  -  CA.product_p(a, A) ; 
assert (all (p ( : )  <  tol)); 

p  =  CA. product (A, a)  -  CA.product_s(A,a)  -  CA.product_p(A,a) ; 
assert (all (p ( : )  <  tol)); 

end 


for  k  =  l:n-l 

f printf  (’ \t°/0d-blade  and  complementary  0/0d-blade\n’  ,k,n-k) ; 

1  construct  complementary  blades 
A  =  e; 

A_inv  =  e; 
for  i  =  1 : k 

a  =  CA. set_vector (v(i , : ) ) ; 

A  =  CA. product (A, a)  ; 

A_inv  =  CA .product (a, A_inv) ; 

end 

B  =  CA . dual (A) ; 

V  =  CA.  product  (A,  B)  ;  %  must  be  a  volume  form  (pseudoscalar) 
7,  check  blade  inversion 

assert (norm (CA. product (A, A_inv)  -  e,inf)  <  tol); 
assert (norm (CA. product (A_inv, A)  -  e,inf)  <  tol); 

1  check  compatibility  with  reversion 

assert (norm(CA. reverse (A)  -  A_inv,inf)  <  tol); 

7.  check  blade  orthogonality 

assert (norm (CA .product_s (A, B) , inf )  <  tol); 

assert (norm (CA .product_p (A, B)  -  V,inf)  <  tol); 

7.  random  vectors 

a  =  CA. set_vector(2.0*rand(N,n)  -  1.0); 


7.  projection  and  rejection 
p  =  CA. product (CA.product_s(a, A) ,A_inv) ; 
r  =  CA. product (CA.product_p(a, A) , A_inv) ; 


7.  check  identities 
assert (norm(p  +  r  -  a, inf)  <  tol); 
assert (norm (CA . product_s (p , r) , inf ) 
assert (norm(CA .product_p(p, A) , inf) 
assert (norm(CA .product_p(r ,B) , inf) 
assert (norm (CA . product_s (p , B) , inf ) 
assert (norm (CA . product_s (r , A) , inf ) 
assert (norm(CA .product_p(a, V) , inf) 

end 


°/0  sum  decomposition 

<  tol)  ;  °/0  orthogonality 

<  tol);  7,  inclusion  of  projected 

<  tol);  /  inclusion  of  rejected 

<  tol);  7,  normality  of  projected 

<  tol) ;  7o  normality  of  rejected 

<  tol)  ;  7«  complete  vector  space 
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end 


%===== - ========== - - - ===== - 

function  test_identity_l (n,tol,N) 
fprintf ( ’Test  of  identity:  a~b  =  -b~a\n’); 

CA  =  cliff ord_algebra(n) ; 

7#  random  vectors 

a  =  CA.set_vector(2.0*rand(N,n)  -  1.0); 
b  =  CA . set_vector (2 . 0*rand(N,n)  -  1.0); 

7,  check  antisymmetry 

q  =  CA. product _p (a, b)  +  CA. product _p(b, a) ; 

assert (all (q( : )  <  tol)); 

end 


%======================================================================== 

function  test_identity_2(n, tol,N) 

fprintf ( ’Test  of  identity:  a. (b~c)  =  (a.b)c  -  (a.c)b\n’); 

CA  =  cliff ord_algebra(n) ; 

7,  random  vectors 

a  =  CA.set_vector(2.0*rand(N,n)  -  1.0); 
b  =  CA . set_vector (2 . 0*rand(N,n)  -  1.0); 
c  =  CA.set_vector(2.0*rand(N,n)  -  1.0); 

%  check  the  identity 

abcl  =  CA.product_s(a,CA.product_p(b,c)) ; 

abc2  =  CA. product (CA.product_s(a,b) ,c)  -  CA. product (CA.product_s(a,c) ,b) ; 

assert (norm (abcl  -  abc2,inf)  <  tol); 

end 


%================ - ====================================== - ====== 

function  test_identity_3(n, tol,N) 

fprintf ( ’Test  of  identity:  (a~b) . (c~d)  =  (a.d)(b.c)  -  (a. c) (b . d) \n’ ) ; 
CA  =  cliff ord_algebra(n) ; 


7*  random  vectors 
a  =  CA.set_vector(2.0*rand(N,n) 
b  =  CA.set_vector(2.0*rand(N,n) 
c  =  CA.set_vector(2.0*rand(N,n) 
d  =  CA.set_vector(2.0*rand(N,n) 


1.0); 

1.0); 

1.0); 

1.0); 


7.  check  the  identity 
ab  =  CA.product_p(a,b) ; 
cd  =  CA.product_p(c,d) ; 
ab_cd  =  CA.product_s(ab,cd) ; 
a_d  =  CA .product_s (a,d) ; 
b_c  =  CA .product_s (b, c) ; 
a_c  =  CA .product_s (a, c) ; 
b_d  =  CA.product_s(b,d) ; 

assert (norm (ab_cd  -  (a_d.*b_c  -  a_c . *b_d) , inf )  <  tol); 
end 


•/,============================================= 

function  test_identity_4(n, tol,N) 

fprintf ( ’Test  of  identity:  a'"C''b  =  -b~C~a\n’); 

CA  =  cliff ord_algebra(n) ; 
d  =  CA. dimension; 

7.  random  vectors 

a  =  CA.set_vector(2.0*rand(N,n)  -  1.0); 
b  =  CA . set_vector (2 . 0*rand(N,n)  -  1.0); 
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7.  random  multivectors 
c  =  2.0*rand(N,d)  -  1.0; 

7#  check  antisymmetry 

acb  =  CA.product_p(a,CA.product_p(c,b) ) ; 
bca  =  CA. product_p(b,CA. product _p(c, a) ) ; 
assert (norm (acb  +  bca, inf)  <  tol) ; 
end 


•/,============================================================ 

function  test_identity_5(n,tol) 

fprintf ( ’Test  of  identity:  a_l~ . . . ~a_n  =  det (a) e_{l . . .n}\n’ ) ; 

CA  =  cliff ord_algebra(n) ; 
d  =  CA. dimension; 

7,  random  basis  vectors 
v  =  2 . 0*rand(n,n)  -  1.0; 

7*  check  the  identity 
p  =  CA . set_scalar (1 . 0) ; 
for  i  =  l:n 

p  =  CA.product_p(p,CA.set_vector(v(i, : ))) ; 

end 

q  =  det (v) ; 

assert (abs(p(d)  -  q)  <  tol); 
end 


%================================================== 

function  test_cross_product (tol ,N) 

fprintf ( ’Test  of  cross  product:  a  x  b  =  *(a~b)\n’); 

CA  =  cliff ord_algebra(3) ; 

7.  random  vectors 
a  =  2 . 0*rand(N,3)  -  1.0; 
b  =  2 . 0*rand(N,3)  -  1.0; 

7.  check  the  identity 

abl  =  CA.set_vector(cross(a,b,2)) ; 

a  =  CA. set_vector (a) ; 

b  =  CA. set_vector (b) ; 

ab2  =  CA. dual (CA. product _p (a, b) ) ; 

assert (norm (abl  -  ab2,inf)  <  tol); 

end 


•/„================================ 

function  test_complex_number (tol ,N) 
fprintf ( ’Complex  number  test\n’); 

CA  =  cliff ord_algebra(3) ; 

7.  random  complex  numbers 

zl  =  2 . 0*rand(N , 2)  -  1.0; 

cl  =  [zl(: ,1) ,zeros(N,6) ,zl(: ,2)] ; 

zl  =  complex(zl(: ,1) ,zl(: ,2)) ; 

z2  =  2 . 0*rand(N , 2)  -  1.0; 

c2  =  [z2( : , 1) ,zeros(N,6) ,z2( : ,2)] ; 

z2  =  complex(z2( : , 1) ,z2( : ,2) ) ; 

7,  check  the  products 
c  =  CA. product (cl, c2) ; 
z  =  zl.*z2; 

cz  =  [real (z) , zeros (N, 6) ,imag(z)] ; 
assert (norm (cz  -  c,inf)  <  tol); 
end 
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•/,====================================== 

function  test_central_subalgebra(tol ,N) 
fprintf ( ’Central  subalgebra  test\n’); 

CA  =  cliff ord_algebra(3) ; 

°/t  random  complex  numbers 
z  =  2.0*rand(N,2)  -  1.0; 
z  =  [z( : , 1) ,zeros(N,6) ,z( : ,2)]  ; 

7.  random  multivectors 
c  =  2 . 0*rand(N,8)  -  1.0; 

7.  check  commutativity 

zc  =  CA. product (z, c) ; 

cz  =  CA. product (c,z) ; 

assert (norm (zc  -  cz,inf)  <  tol) ; 

end 


•/,===================================== 

function  test_bicomplex_number (tol ,N) 
fprintf ( ’Bicomplex  number  test\n’); 

CA  =  cliff ord_algebra(3) ; 

7.  random  bi complex  numbers 
a  =  2 . 0*rand(N,4)  -  1.0; 
a  =  [a( : , 1 :2) ,zeros(N,4) ,a( : ,3:4)] ; 
b  =  2 . 0*rand(N,4)  -  1.0; 
b  =  [b(: ,1:2) ,zeros(N,4) ,b(: ,3:4)]  ; 

7.  check  commutativity 
q  =  CA. product (a, b)  -  CA. product (b, a) ; 
assert (all (q( : )  <  tol)); 
end 


%================ - ============================== 

function  test_quaternion(tol ,N) 
fprintf (’ Quaternion  test\n’); 

CA  =  cliff ord_algebra(3) ; 
e  =  CA. set_scalar (ones (N,  1) ) ;  7«  units 

7,  random  quaternions  and  their  inverses 
q  =  CA.even(2.0*rand(N,8)  -  1.0); 
q_inv  =  CA. reverse (q) ; 
s  =  CA. product _s(q,q_inv) ; 
s  (  :  ,  1)  =  1.0./s(:,l); 
q_inv  =  CA. product (s, q_inv) ; 

7.  check  quaternion  inversion 

assert (norm (CA. product (q,q_inv)  -  e,inf)  <  tol); 
assert (norm (CA. product (q_inv,q)  -  e,inf)  <  tol); 

7,  random  vectors  and  their  vector/scalar  products 
a  =  CA.set_vector(2.0*rand(N,3)  -  1.0); 
b  =  CA.set_vector(2.0*rand(N,3)  -  1.0); 
c  =  CA.dual(CA.product_p(a,b)) ; 
s  =  CA .product_s (a,b) ; 

7,  rotated  vectors  and  their  vector/scalar  products 
al  =  CA. product (q,CA. product (a, q_inv)) ; 
bl  =  CA. product (q,CA. product (b,q_inv)) ; 
cl  =  CA.dual(CA.product_p(al,bl)) ; 
si  =  CA.product_s(al ,bl) ; 
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7,  check  space  invariance 

a2  =  CA.set_vector(CA.get_vector(al)) ; 

assert (norm (a2  -  al,inf)  <  tol) ; 

b2  =  CA.set_vector(CA.get_vector(bl)) ; 

assert (norm (b2  -  bl,inf)  <  tol); 

%  check  isometric  property 
assert (norm(sl  -  s,inf)  <  tol); 

7,  check  orientation  preserving  property 
c2  =  CA. product (q,CA. product (c,q_inv)) ; 
assert (norm (c2  -  cl, inf)  <  tol); 
end 


•/,=========================================================== 

function  test_perf ormance (N) 
fprintf ( ’Performance  test\n’); 

CA  =  cliff ord_algebra(4) ; 

CA. display () ; 
d  =  CA. dimension; 

7,  random  multivectors 
A  =  2.0*rand(N,d)  -  1.0; 
b  =  2.0*rand(l,d)  -  1.0; 

fprintf ( ’Array  size:  7od\n’,N); 

7.  vectorized  multiplication 
tic ; 

CA . product ( A , b) ; 
vector_time  =  toe; 

7.  serial  multiplication 
tic ; 

for  i  =  1:N 

a  =  A(i , : ) ; 

CA. product (a, b) ; 

end 

serial_time  =  toe; 

fprintf  ( ’Performance  factor:  7»g\n’ , serial_time/vector_time)  ; 
end 
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Appendix  D  Listing  of  geometric  algebra  test  .m 

‘/.GEOMETRIC.  ALGEBRA_TEST  (N ,  tol) 

7. 

7. Unit  tests  of  application  of  geometric  algebra  to  projective  geometry. 

7. 

/.Copyright  (C)  2014  Defence  Science  and  Technology  Organisation 

7. 

/.Created  by  Leonid  K.  Antanovskii 
function  geometric_algebra_test (N , tol) 
fprintf ( ’Testing  geometric  algebra\n’); 
if  nargin  <  2 

rng( ’default ’ ) ; 
tol  =  1.0e-13; 
if  nargin  <  1 
N  =  100000; 

end 

end 

test_projective_space2D(tol ,N) ; 
test_projective_space3D(tol ,N) ; 
test_pappus_theorem(tol ,N) ; 
test_desargues_theorem(tol ,N) ; 
test_plucker_coordinates(tol ,N) ; 
test_epipolar_geometry (tol ,N) ; 
test_fundamental_map(tol ,N) ; 
test_point_reconstruction(tol ,N) ; 
test_line_reconstruction(tol ,N) ; 

end 


•/,=============================================== 

function  test_projective_space2D(tol ,N) 
fprintf (’ Incidence  test  in  projective  plane\n’); 

CA  =  cliff ord_algebra(3) ; 

7.  original  points 

XI  =  CA.set_vector(2.0*rand(N,3)  -  1.0); 

X2  =  CA.set_vector(2.0*rand(N,3)  -  1.0); 

X3  =  CA.set_vector(2.0*rand(N,3)  -  1.0); 

/.  lines  passing  through  point  pairs 
L12  =  CA. product _p (XI, X2) ; 

L13  =  CA. product _p (XI, X3) ; 

L23  =  CA.product_p(X2,X3) ; 

7.  intersection  points 
Y1  =  CA.product.r (L12,L13) ; 

Y2  =  CA . product.r (L12 , L23) ; 

Y3  =  CA . product.r (L13 , L23) ; 

7#  comparison 
q  =  CA .product_p(Xl , Yl) ; 
assert (all (q( : )  <  tol)); 
q  =  CA.product_p(X2,Y2) ; 
assert (all (q( : )  <  tol)); 
q  =  CA.product_p(X3,Y3) ; 
assert (all (q( : )  <  tol)); 
end 


%============= - ========== - =========== 

function  test_projective_space3D(tol ,N) 
fprintf (’ Incidence  test  in  projective  space\n’); 

CA  =  cliff ord_algebra(4) ; 
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7#  original  points 

XI  =  CA.set_vector(2.0*rand(N,4)  -  1.0); 

X2  =  CA.set_vector(2.0*rand(N,4)  -  1.0); 

X3  =  CA.set_vector(2.0*rand(N,4)  -  1.0); 

X4  =  CA.set_vector(2.0*rand(N,4)  -  1.0); 

7.  planes  passing  through  point  triplets 
P123  =  CA.product_p(CA.product_p(Xl,X2) ,X3) ; 
P124  =  CA.product_p(CA.product_p(Xl,X2) ,X4) ; 
P134  =  CA.product_p(CA.product_p(Xl,X3) ,X4) ; 
P234  =  CA.product_p(CA.product_p(X2,X3) ,X4) ; 

7*  intersection  points 

Y1  =  CA.product_r (CA.product_r (P123,P124) ,P134) 
Y2  =  CA . product_r (CA . product_r (P123 , P124) , P234) 
Y3  =  CA.product_r (CA.product_r (P123,P134) ,P234) 
Y4  =  CA.product_r (CA.product_r (P124,P134) ,P234) 

7.  comparison 
q  =  CA .product _p (XI ,Y1) ; 
assert (all (q( : )  <  tol)); 
q  =  CA.product_p(X2,Y2) ; 
assert (all (q( : )  <  tol)); 
q  =  CA.product_p(X3,Y3) ; 
assert (all (q( : )  <  tol)); 
q  =  CA.product_p(X4,Y4) ; 
assert (all (q( : )  <  tol)); 
end 


%======================================================== 

function  test_pappus_theorem(tol ,N) 
fprintf ( ’Pappus ’  3  theorem  test\n’); 

CA  =  cliff ord_algebra(3) ; 

7,  first  triplet  of  collinear  points 
A1  =  CA.set_vector(2.0*rand(N,3)  -  1.0); 

B1  =  CA.set_vector(2.0*rand(N,3)  -  1.0); 

Cl  =  rand*Al  +  rand*Bl; 

7.  check  collinearity 

q  =  CA.product_p(CA.product_p(Al ,B1) ,C1) ; 
assert (all (q( : )  <  tol)); 

7.  second  triplet  of  collinear  points 
A2  =  CA.set_vector(2.0*rand(N,3)  -  1.0); 

B2  =  CA.set_vector(2.0*rand(N,3)  -  1.0); 

C2  =  rand*A2  +  rand*B2; 

7,  check  collinearity 

q  =  CA . product_p (CA . product_p (A2 , B2) , C2) ; 
assert (all (q( : )  <  tol)); 

7.  intersection  points 

A  =  CA.product_r(CA.product_p(Bl,C2) ,CA.product_p(B2,Cl)) 
B  =  CA.product_r(CA.product_p(Cl,A2) ,CA.product_p(C2,Al)) 
C  =  CA.product_r(CA.product_p(Al,B2) ,CA.product_p(A2,Bl)) 

7.  check  collinearity  of  the  intersection  points 
q  =  CA.product_p(CA.product_p(A,B) ,C) ; 
assert (all (q( : )  <  tol)); 
end 


•/,===================================== 

function  test_desargues_theorem(tol ,N) 
fprintf ( 3 Desargues 3  3  theorem  test\n ’ ) ; 


UNCLASSIFIED 


43 


DSTO-TR  3021 


UNCLASSIFIED 


CA  =  cliff ord_algebra(3) ; 

7,  triangle  vertices  in  projective  plane 
A1  =  CA.set_vector(2.0*rand(N,3)  -  1.0); 

B1  =  CA.set_vector(2.0*rand(N,3)  -  1.0); 

Cl  =  CA.set_vector(2.0*rand(N,3)  -  1.0); 

A2  =  CA.set_vector(2.0*rand(N,3)  -  1.0); 

B2  =  CA.set_vector(2.0*rand(N,3)  -  1.0); 

C2  =  CA.set_vector(2.0*rand(N,3)  -  1.0); 

7.  intersection  points  of  triangle  sides 

A  =  CA.product_r(CA.product_p(Bl,Cl) ,CA.product_p(B2,C2)) 
B  =  CA.product_r(CA.product_p(Cl,Al) ,CA.product_p(C2,A2)) 
C  =  CA.product_r(CA.product_p(Al,Bl) ,CA.product_p(A2,B2)) 

7.  test  for  collinearity  of  points 
p  =  CA . product_p (CA . product_p (A , B) , C) ; 

7,  lines  through  triangle  vertices 
La  =  CA.product_p(Al ,A2) ; 

Lb  =  CA.product_p(Bl ,B2) ; 

Lc  =  CA.product_p(Cl ,C2) ; 

7.  test  for  concurrency  of  lines 
q  =  CA . product_r (CA . product_r (La , Lb) , Lc) ; 

7.  non-degenerate  transformations 
p  =  CA.dual(p); 

11  =  CA.product_p(CA.product_p(Al ,B1) ,C1) ; 

12  =  CA . product_p (CA . product_p (A2 , B2) , C2) ; 
s  =  -CA. product (il, i2) ; 

q  =  CA .product (s ,q) ; 

7,  (p  =  0  <=>  q  =  0)  follows  from  p  =  q 

assert (norm (p  -  q,inf)  <  tol) ; 

end 


•/,======================================= 

function  test_plucker_coordinates (tol , N) 
fprintf ( ’Plucker  coordinates  test\n’); 

CA  =  cliff ord_algebra(4) ; 

7.  random  points 

XI  =  CA.set_vector(2.0*rand(N,4)  -  1.0); 
X2  =  CA.set_vector(2.0*rand(N,4)  -  1.0); 

7.  lines  through  points 
L  =  CA.product_p(Xl ,X2) ; 

7.  check  Klein  quadric 
q  =  klein_quadric(L) ; 
assert(all(q  <  tol)); 

7,  random  planes 
PI  =  CA. dual (XI); 

P2  =  CA . dual (X2) ; 

7,  lines  through  planes 
L  =  CA.product_r (PI ,P2) ; 

7,  check  Klein  quadric 
q  =  klein_quadric(L) ; 
assert(all(q  <  tol)); 
end 
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•/,================================================ 

function  test_epipolar_geometry (tol ,N) 
fprintf ( ’Epipolar  geometry  test\n’); 

CA  =  cliff ord_algebra(4) ; 

7,  first  camera  centre  and  projection  plane 
Cl  =  CA.set_vector(2.0*rand(l,4)  -  1.0); 

PI  =  CA.dual(CA.set_vector(2.0*rand(l,4)  -  1.0)); 

7.  second  camera  centre  and  projection  plane 
C2  =  CA.set_vector(2.0*rand(l,4)  -  1.0); 

P2  =  CA.dual(CA.set_vector(2.0*rand(l,4)  -  1.0)); 

7,  base  line  and  epipoles 
B  =  CA.product_p(Cl ,C2) ; 

01  =  CA.product_r (B,P1) ; 

02  =  CA.product_r (B,P2) ; 

7,  random  points 

X  =  CA . set_vector (2 . 0*rand(N,4)  -  1.0); 

7,  planes  through  the  baseline  and  the  points 
P  =  CA.product_p(B,X) ; 

7.  projection  rays 

LI  =  CA.product_p(Cl ,X) ; 

L2  =  CA.product_p(C2,X) ; 

7,  projection  points 

XI  =  CA.product_r (PI ,L1) ; 

X2  =  CA.product_r(P2,L2) ; 

7,  epipolar  lines 

El  =  CA.product_r (P,P1) ; 

E2  =  CA.product_r (P,P2) ; 

7.  check  incidence 
q  =  CA.product_p(El ,X1) ; 
assert (all (q( : )  <  tol)); 
q  =  CA.product_p(E2,X2) ; 
assert (all (q( : )  <  tol)); 
q  =  CA.product_p(El ,01) ; 
assert (all (q( : )  <  tol)); 
q  =  CA.product_p(E2,02) ; 
assert (all (q( : )  <  tol)); 
end 


%================================================ 

function  test_f undamental_map (tol , N) 
fprintf ( ’Fundamental  map  test\n’); 

CA  =  cliff ord_algebra(4) ; 

7,  first  camera  centre  and  projection  plane 
Cl  =  CA.set_vector(2.0*rand(l,4)  -  1.0); 

PI  =  CA.dual(CA.set_vector(2.0*rand(l,4)  -  1.0)); 

7,  second  camera  centre  and  projection  plane 
C2  =  CA.set_vector(2.0*rand(l,4)  -  1.0); 

P2  =  CA.dual(CA.set_vector(2.0*rand(l,4)  -  1.0)); 

7,  base  line  and  epipoles 
B  =  CA.product_p(Cl ,C2) ; 

01  =  CA.product_r (B,P1) ; 

02  =  CA.product_r (B,P2) ; 
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7,  planes  through  the  baseline  and  random  points 
P  =  CA.product_p(B,CA. set_vector(2.0*rand(N,4)  -  1.0)); 

%  epipolar  lines  induced  by  the  planes 
El  =  CA.product_r (P,P1) ; 

E2  =  CA.product_r (P,P2) ; 

7*  check  incidence  with  the  epipoles 
assert (norm (CA.product_p (01, El) , inf)  <  tol) ; 
assert (norm (CA.product_p (02, E2) , inf )  <  tol); 

7,  action  of  the  fundamental  map 

FI  =  CA.product_r(CA.product_p(C2,E2) ,P1) ;  %  E2  ->  FI 
F2  =  CA.product_r (CA.product_p(Cl ,E1) ,P2) ;  %  El  ->  F2 

7.  check  incidence  with  the  epipoles 
assert (norm (CA.product_p (01, FI) , inf)  <  tol); 
assert (norm (CA.product_p (02, F2) , inf )  <  tol); 

7.  check  line  equalities 
for  i  =  1:N 
A  =  [ 

line_matrix(El (i , : ) ) 
line_matrix(Fl (i , : ) ) 

]; 

assert (rank(A, tol)  ==  2); 

end 

for  i  =  1:N 
A  =  [ 

line_matrix(E2(i , :)) 
line_matrix(F2(i , :)) 

]; 

assert (rank(A, tol)  ==  2); 

end 

end 

%===== - ========== - - - ======== - =========== 

function  test_point_reconstruction(tol ,N) 
fprintf ( ’Point  reconstruction  test\n’); 

CA  =  cliff ord_algebra(4) ; 

7,  first  camera  centre  and  projection  plane 
Cl  =  CA.set_vector(2.0*rand(l,4)  -  1.0); 

PI  =  CA.dual(CA.set_vector(2.0*rand(l,4)  -  1.0)); 

7.  second  camera  centre  and  projection  plane 
C2  =  CA.set_vector(2.0*rand(l,4)  -  1.0); 

P2  =  CA.dual(CA.set_vector(2.0*rand(l,4)  -  1.0)); 

7.  random  points  to  project  and  recover 

X  =  CA.set_vector(2.0*rand(N,4)  -  1.0); 

7,  projection  rays 

LI  =  CA.product_p(Cl ,X) ; 

L2  =  CA.product_p(C2,X) ; 

7,  projection  points 

XI  =  CA.product_r (PI ,L1) ; 

X2  =  CA.product_r(P2,L2) ; 

7,  check  inclusion  using  both  outer  products 
assert (norm(CA.product_p(Xl , PI) ,inf)  <  tol); 
assert (norm (CA. product_p (X2, P2) , inf )  <  tol); 
assert (norm(CA.product_r(Xl , PI) ,inf)  <  tol); 
assert (norm (CA. product_r (X2, P2) , inf )  <  tol); 
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7,  reconstructed  projection  rays 
LI  =  CA.product_p(Cl ,X1) ; 

L2  =  CA.product_p(C2,X2) ; 

7*  check  signed  crossing  values 
q  =  CA.product_p(Ll ,L2) ; 
assert (all (q( : )  <  tol)); 

7,  check  incidence  with  ray  1 
q  =  CA.product_p(Ll ,X) ; 
assert (all (q( : )  <  tol)); 

7,  check  incidence  with  ray  2 
q  =  CA.product_p(L2,X) ; 
assert (all (q( : )  <  tol)); 

7.  line  intersection  points 
v  =  zeros (N,4) ; 
for  i  =  1:N 
A  =  [ 

line_matrix(Ll (i , : ) ) 
line_matrix(L2(i , :)) 

]; 

[~,D,V]  =  svd(A,0) ; 
assert (abs(D(4, 4) )  <  tol); 
v(i , : )  =  V( : ,4) ’ ; 

end 

Y  =  CA. set_vector (v) ; 

7.  check  recovered  triangulation  points 
q  =  CA.product_p(X,Y) ; 
assert (all (q( : )  <  tol)); 
end 


%========================================================== 

function  test_line_reconstruction(tol ,N) 
fprintf ( ’Line  reconstruction  test\n’); 

CA  =  cliff ord_algebra(4) ; 

7,  first  camera  centre  and  projection  plane 
Cl  =  CA.set_vector(2.0*rand(l,4)  -  1.0); 

PI  =  CA.dual(CA.set_vector(2.0*rand(l,4)  -  1.0)); 

7,  second  camera  centre  and  projection  plane 
C2  =  CA.set_vector(2.0*rand(l,4)  -  1.0); 

P2  =  CA.dual(CA.set_vector(2.0*rand(l,4)  -  1.0)); 

7,  random  lines  to  project  and  recover 
L  =  CA.product_p( . . . 

CA.set_vector(2.0*rand(N,4)  -  1.0),... 
CA.set_vector(2.0*rand(N,4)  -  1.0)); 

7,  projected  lines 

LI  =  CA.product_r (PI ,CA.product_p(Cl ,L) ) ; 

L2  =  CA.product_r(P2,CA.product_p(C2,L)) ; 

7.  check  inclusion 

assert (norm (CA.product_r (LI, PI) , inf)  <  tol); 
assert (norm (CA. product_r (L2, P2) , inf )  <  tol); 

7,  reconstructed  lines 

L0  =  CA. product _r(CA. product _p (Cl , LI) ,CA. product _p(C2,L2) ) ; 

7,  check  line  equalities 
for  i  =  1:N 
A  =  [ 
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line_matrix (L ( i , : ) ) 
line_matrix(LO (i , : ) ) 
]; 

assert(rank(A,tol)  ==  2); 

end 

end 


•/,================================== 

function  q  =  klein_quadric(L) 

112  =  L( : ,6) ; 

113  =  L( : ,7); 

114  =  L( : ,8) ; 

123  =  L( : ,9) ; 

124  =  L(: ,10); 

134  =  L(: ,11); 

q  =  112. *134  -  113. *124  +  114. *123; 
end 


function  A  =  line_matrix(L) 
A  =  [ 

L(9) ,-L(7) ,L(6) ,0.0 
L(10) ,-L(8) ,0.0, L(6) 
L(ll) ,0.0, -L(8) ,L(7) 
0.0,L(11) ,-L(10) ,L(9) 

]; 

end 
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